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1.  Summary 


Our  efforts  during  the  project  were  divided  among  two  components  of  study:  ( 1 )  regional 
waveform  modeling  and  hypocenter  determination  and  (2)  joint  satellite  image  and  teleseismic 
location  analysis.  Results  from  each  component  of  our  research  are  summarized  below. 

Three  component  seismic  data  from  a  set  of  presumed  explosions  recorded  by  stations  at 
Bayanaul  and  Karkaralinsk  in  Kazakhstan  were  analyzed  in  order  to  model  the  crustal  structure  of 
the  region  and  to  examine  the  use  of  the  arrival  times  of  secondary  P  phases  in  regional  event 
location.  Data  from  the  first  5  to  10  seconds  of  13  presumed  explosions  were  modeled  with  the 
reflectivity  method.  A  good  fit  to  the  arrival  times  and  amplitudes  in  the  first  five  seconds  of  the  P 
wave  was  obtained  in  the  epicentral  distance  range  of  100  to  300  km.  The  crustal  P  wave  velocity 
model  we  derived  has  an  upper  crustal  velocity  increasing  fairly  rapidly  from  4.5  km/sec  near  the 
surface  to  6.5  km/sec  at  15  km  depth,  then  increasing  more  slowly  to  7.05  km/sec  at  50  km  depth. 
We  used  the  derived  crustal  model  and  the  primary  and  secondary  P  wave  arrival  times  to  relocate 
events  in  the  Kazakhstan  region.  Inclusion  of  the  phase  PmP  substantially  decreases  the  focal 
depth  uncertainty  for  many  events.  All  but  o*e  of  the  events  analyzed  are  concluded  to  be  surface 
explosions;  the  identity  of  the  remaining  event  is  uncertain. 

Using  this  velocity  model  for  the  source  region,  we  extended  our  analysis  to  greater 
distances.  Waveforms  from  nuclear  explosions  at  the  Kazakhstan  test  site  in  the  years  1971-1989 
recorded  at  four  stations  in  the  distance  range  600-1000  km  were  modeled  using  a  modified 
version  of  the  generalized  ray  Cagniard  code  that  includes  separate  source  and  receiver  velocity 
structures  to  account  for  the  effect  of  lateral  heterogeneity.  We  obtained  separate  receiver  models 
for  each  station  region  with  the  source  region  model  obtained  above.  We  match  the  observed  first 
arrival  within  about  0.3  seconds  by  varying  the  crust  and  upper  mantle  structure  at  the  receiver. 
The  secondary  phases  from  the  upper  mantle  are  matched  within  about  1  or  2  seconds  in  the  first 
10  seconds  after  the  first  arrival,  and  we  also  obtain  good  fits  to  the  overall  waveform  envelopes. 

Locations  of  20  nuclear  explosions  from  1987  to  1989  at  the  Balapan  test  site  in 
Kazakhstan  were  derived  with  a  precision  of  about  100  m  from  a  combination  of  time-sequence 
satellite  images  and  teleseismic  epicenter  estimates.  Ground  control  points  for  satellite  image 
rectification  were  obtained  from  information  on  Balapan  explosions  published  by  Bocharov  et  al. 
(1989).  Fresh  disturbances  mapped  in  rectified  SPOT  satellite  images  have  been  associated  with 
individual  explosions.  The  seismically-determined  locations  are  associated  with  formal  error 
estimates  (95%  confidence  ellipses)  that  are  significantly  too  small.  In  15  out  of  19  cases,  the 
events  are  found  from  satellite  imagery  to  lie  outside  these  ellipses.  Our  effort  was  then  extended 
to  the  analysis  of  older  LANDSAT  images.  We  have  been  successful  in  identifying  nearly  all 
events  from  1973  through  1985. 
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2.  Accomplishments 


2.1.  Task  Objectives 

The  goal  of  this  project  has  been  to  address  some  of  the  difficulties  faced  in  regional  and 
teleseismic  event  location  for  seismic  verification  puiposes.  The  data  we  have  analyzed  include 
near-regional  and  far-regional  seismograms,  teleseismic  arrival  times,  and  satellite  images. 

For  regional  event  location,  we  have  been  concerned  with  the  effect  of  source  distance, 
source  depth,  and  regional  velocity  structure  on  the  initial  P  waveform  (about  10  seconds 
following  first-P).  At  near-regional  distances  (less  than  about  500  km),  a  single  plane-layered 
velocity  structure  tends  to  be  adequate  for  waveform  modeling.  In  this  distance  range,  distinct 
secondary  P  phases  often  can  be  observed  and  identified,  and  hence  can  provide  valuable 
constraint  on  source  location  and  depth.  The  data  we  have  analyzed  are  from  the  former 
NRDC/Soviet  Academy  of  Sciences  network  in  Kazakhstan.  Our  objective  has  been  to  confirm  the 
identity  of  secondary  phases  via  waveform  modeling,  improving  the  model  of  velocity  structure  in 
the  vicinity  of  the  Kazakhstan  Test  Site  (KTS)  in  the  process,  and  evaluating  the  contribution  of 
secondary  phases  to  event  location  capability  (see  section  2.4.1).  At  mid-  to  far-regional  distances 
(500  to  2000  km),  a  single  plane-layered  velocity  structure  is  generally  no  longer  adequate  for 
waveform  modeling,  and  secondary  P  phases  are  more  difficult  to  identify.  Our  objectives  have 
been  to  adapt  the  generalized  ray  seismogram  synthesis  technique  to  incorporate  separate  source 
and  receiver  structures,  and  use  this  method  to  model  regional  P  waveforms  and  velocity  structure 
in  Central  Asia,  using  sources  from  KTS  and  vicinity  (see  section  2.4.2).  Using  this  modeling 
capability,  we  begin  to  address  the  question  of  depth  determination  for  seismic  events  recorded  at 
moderate  regional  distances  based  on  initial  P  waveform  characteristics. 

For  teleseismic  location,  our  basic  objective  is  to  establish  a  framework  for  evaluation  of 
teleseismic  location  capability  by  deriving  "ground  truth"  information  for  the  majority  of  nuclear 
explosions  with  mb  >  5  at  the  Balapan  (Shagan  River)  area  of  KTS.  Initially,  our  efforts  involved 
the  use  of  SPOT  images  for  recent  (1987  to  1989)  nuclear  explosions  Firm  identification  was 
made  of  the  sites  of  18  of  the  20  events  from  this  period,  plus  provisional  identification  of  the 
remaining  2  (see  section  2.4.3).  We  then  turned  our  attention  to  LANDSAT  images  covering  the 
time  period  1974  to  1982.  Despite  the  lower  resolution  of  the  LANDSAT  data  (80  m,  versus  10  or 
20  m  for  SPOT),  we  again  succeeded  in  firmly  identifying  the  vast  majority  of  the  sites  of  nuclear 
explosions  in  this  period,  and  using  the  1982  LANDSAT  and  1986  SPOT  images,  nearly  all  of  the 
event  sites  between  1982  and  1986  (see  section  2.4.4).  With  this  identification  effort  completed  to 
the  extent  possible,  these  "ground  truth"  results  can  be  put  to  future  use  for  a  variety  of  purposes. 
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2.2.  Technical  Problem 

Seismic  event  location  remains  a  fundamental  component  of  monitoring  efforts  related  to 
verifying  nuclear  test  ban  or  test  limitation  treaties.  Event  location,  especially  depth  determination, 
is  important  both  as  a  basis  for  discrimination  and  as  a  starting  point  for  the  analysis  of  wave 
propagation  and  attenuation.  The  earthquake  location  problem  is  relatively  well  understood  on  a 
theoretical  basis.  However,  it  can  be  expected  that  event  location  will  be  a  non-trivial  problem  for 
in-country  regional  networks  and  global  non-proliferation  monitoring.  The  following  paragraphs 
are  adapted  from  the  contribution  of  the  PI  to  DARPA's  Seismic  Identification  Workshop  report. 

Considerable  progress  has  been  achieved  on  the  problem  of  epicenter  determination  for 
small  events  at  regional  distances  with  sparse  networks  and/or  arrays.  One  of  the  keys  to 
improved  results  is  the  incorporation  of  arrival  azimuth  and  slowness  information  in  the  location 
procedure.  A  number  of  recent  studies  have  demonstrated  the  capability  for  epicenter  accuracy  of  a 
few  tens  of  kilometers  at  regional  distances.  It  is  also  apparent  that  location  accuracy  can  be 
improved  if  information  is  available  on  crust  and  mantle  structure  between  the  sources  and 
receivers  ,  as  well  as  arrival  azimuth  anomalies  at  the  receivers.  Thus  there  is  a  need  to  develop 
and  apply  methods  for  calibrating  velocity  structure  in  (possibly  remote)  source  regions  and 
evaluating  near-receiver  wave  propagation  anomalies.  Travel  time  and  azimuth  uncertainties  must 
be  reduced  substantially  in  order  to  achieve  a  level  of  location  accuracy  (10  to  20  km)  that  would 
permit  identification  of  potential  surface  sources  of  seismic  events  in  satellite  imagery. 

The  problem  of  depth  determination  has  proven  to  be  much  more  difficult  than  epicenter 
determination.  At  short  distances  (<  500  km),  Rg  excitation  is  often  diagnostic  of  near-surface 
sources,  but  Rg  does  not  usually  propagate  to  large  distance.  At  near-regional  distances  (500  to 
1000  km),  empirical  and  theoretical  works  indicate  that  reasonable  depth  constraint  (±  5  km)  can  be 
obtained  if  arrival  times  of  secondary  phases  (Pg,  Sn,  Lg)  are  available  and  lithospheric  structure  is 
accurately  known.  However,  at  far-regional  distances  (1000  to  2000  km),  the  combination  of 
lateral  heterogeneity,  crustal  attenuation,  and  crustal  reverberation  is  likely  to  prevent  routine  use  of 
secondary  phase  times  other  than  Lg,  which  has  the  greatest  timing  (and  modeling)  uncertainty. 
Instead,  modeling  of  body-wave  depth  phases  (pP,  sP,  etc.)  would  appear  to  be  necessary.  Depth 
phases  are  often  detected  reliably  by  arrays,  but  methods  to  identify  them  from  sets  of  single¬ 
station  recordings  are  needed.  The  goal  in  depth  estimation  should  be  an  ability  to  determine  with 
90%  certainty  whether  or  not  a  particular  event  occurred  within  1  kilometer  of  the  surface. 

Progress  in  epicenter  and  depth  determination  will  require  both  theoretical  and  empirical 
investigations,  the  latter  based  on  data  for  which  "ground  truth"  is  available.  The  limitations  of 
existing  and  new  approaches  need  to  be  determined  for  appropriate  nonproliferation  scenarios.  In 
particular,  it  is  critical  to  know  how  much  additional  information  might  be  required  for  successful 
location  and  depth  determination  of  suspicious  events  in  new  (previously  uncalibrated)  regions. 
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2.3.  General  Methodology 


Our  methodology  for  this  project  involves  3  primary  elements:  (1)  regional  seismic 
waveform  modeling;  (2)  teleseismic  joint  epicenter  determination,  and  (3)  satellite  image  analysis. 

Three  component  seismic  data  from  a  set  of  presumed  explosions  recorded  by  stations  at 
Bayanaul  and  Karkaraiinsk  in  Kazakhstan  were  analyzed  in  order  to  model  the  crustal  structure  of 
the  region  and  to  examine  the  use  of  the  arrival  times  of  secondary  P  phases  in  regional  event 
location.  Data  from  the  first  5  to  10  seconds  of  13  presumed  explosions  were  modeled  with  the 
reflectivity  method,  and  a  new  regional  velocity  structure  was  obtained.  Using  this  velocity  model 
for  the  source  region,  we  extended  our  analysis  to  stations  at  greater  distances.  Waveforms  from 
nuclear  explosions  at  the  Kazakhstan  test  site  in  the  years  197 1  to  1989  recorded  at  four  stations  in 
the  distance  range  600  to  1000  km  were  modeled  using  a  modified  version  of  the  generalized  ray 
Cagniard  code  that  includes  separate  source  and  receiver  velocity  structures  to  account  for  the  effect 
of  lateral  heterogeneity.  We  obtained  separate  receiver  models  for  each  station  region  with  the 
source  region  model  obtained  above. 

For  our  teleseismic  location  efforts,  we  utilized  the  algorithms  JED  and  MLOC  for  joint 
epicenter  determination.  We  examined  the  sensitivity  of  the  locations  to  the  use  of  single  or 
multiple  master  events.  Our  goal  was  to  establish  an  estimate  of  the  level  of  uncertainty  that  might 
be  associated  with  these  location  estimates,  so  that  this  uncertainty  would  be  taken  into  account 
when  time- sequence  satellite  images  were  subsequently  analyzed. 

Ground  control  points  for  SPOT  satellite  image  rectification  were  obtained  from 
information  on  Balapan  explosions  published  by  Bocharov  et  al.  (1989).  Locations  of  20  nuclear 
explosions  from  1987  to  1989  at  the  Balapan  test  site  in  Kazakhstan  were  derived  with  a  precision 
of  about  100  m  from  a  combination  of  time-sequence  SPOT  satellite  images  and  teleseismic 
epicenter  estimates.  Our  effort  was  then  extended  to  the  analysis  of  older  LANDSAT  images. 
Enhancement  techniques  were  evaluated  in  order  to  maximize  the  visibility  of  shot  point  features. 

2.4.  Technical  Results 

The  project  duration  was  approximately  28  months,  8/21/90  to  12/31/92.  Our  efforts 
during  the  period  were  divided  among  two  components  of  study:  (1)  regional  waveform  modeling 
and  hypocenter  determination  and  (2)  joint  satellite  image  and  teleseismic  location  analysis. 

Progress  achieved  in  these  two  project  components  is  described  in  the  following  four  subsections. 

2.4.1.  Seismic  velocity  structure  and  event  relocation 
in  Kazakhstan  from  secondary  P  phases 
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SEISMIC  VELOCITY  STRUCTURE  AND  EVENT  RELOCATION 
IN  KAZAKHSTAN  FROM  SECONDARY  P  PHASES 

By  H.  R.  Quin  and  C.  H.  Thurber 
Abstract 

Three-component  seismic  data  from  a  set  of  presumed  explosions  recorded 
by  stations  at  Bayanaul  and  Karkaralinsk  in  Kazakhstan  were  analyzed  in  order 
to  model  the  crustal  structure  of  the  region  and  to  examine  the  use  of  the 
arrival  times  of  secondary  P  phases,  primarily  PmP,  in  regional  event  location. 
Polarization  analysis  aided  in  the  identification  of  the  secondary  phases. 
Low-pass  filtered  data  (4-Hz  corner)  from  the  first  5  to  10  sec  of  13  presumed 
explosions  were  modeled  with  the  reflectivity  method.  The  two  chemical 
explosions  in  1987  provided  a  check  on  accuracy,  as  their  locations  and  origin 
times  are  accurately  known.  A  good  fit  to  the  arrival  times  and  amplitudes  in 
the  first  5  sec  of  the  P  wave  (Pn,  Pg,  and  PmP)  was  obtained  in  the  epicentral 
distance  range  of  100  to  300  km.  Beyond  300  km,  the  simple  layered  model 
was  not  adequate  to  model  the  PmP  arrival. 

The  crustal  P- wave  velocity  model  we  derived  has  ah  upper  crustal  velocity 
increasing  fairly  rapidly  from  4.5  km  /  sec  near  the  surface  to  6.5  km  /  sec  at 
15-km  depth,  then  increasing  more  slowly  to  7.05  km  /  sec  at  50-km  depth.  The 
observed  difference  in  the  arrival  times  of  the  phases  Pg,  PmP,  and  Pn  in  the 
range  between  100-  and  250-km  distance  required  a  relatively  sharp  transition 
at  the  crust  mantle  boundary.  The  model  is  generally  similar  to  previous 
estimates  of  P  velocity  structure  In  the  region,  though  with  a  gentler  gradient 
in  the  upper  crust  and  a  steeper  gradient  In  the  lower  crust.  We  used  the 
derived  crustal  model  and  the  primary  and  secondary  P-wave  arrival  times  to 
relocate  events  in  the  Kazakhstan  region.  Inclusion  of  the  phase  PmP  substan¬ 
tially  decreases  the  focal  depth  uncertainty  for  many  of  the  events.  All  but  one 
of  the  events  analyzed  are  concluded  to  be  surface  explosions;  the  identity  of 
the  remaining  event  is  uncertain. 

Introduction 

Seismic  event  location  at  regional  distances  has  traditionally  relied  on  P- wave 
first  arrivals  at  a  large  number  of  stations.  However,  in  the  case  of  a  sparse 
regional  network  or  regional  arrays,  event  location  often  must  be  done  by 
incorporating  seismic  wave  azimuths  and/or  timing  of  both  first  and  secondary 
arrivals.  Previous  work  on  near-regional  data  from  Kazakhstan  (Thurber  et  al., 
1989;  Li  and  Thurber,  1991)  indicates  that  secondary  P  phases  can  be  observed 
in  many  seismograms.  Secondary  P  phases  have  been  modeled  in  many  studies 
in  other  areas  (e.g.,  Helmberger  and  Engen,  1980;  Langston,  1982;  Holt  and 
Wallace,  1989;  Vogfjord  and  Langston,  1991),  especially  long-period  Pnl,  but  to 
our  knowledge  these  phases  have  not  been  used  for  regional  event  location 
studies  except  for  Thurber  et  al.  (1989)  and  Li  and  Thurber  (1991).  We  have 
undertaken  a  study  to  determine  whether  or  not  secondary  phases  can  be 
accurately  identified,  using  a  combination  of  polarization  analysis  and  seismo¬ 
gram  synthesis,  and  then  utilized  to  improve  event  locations. 

This  study  is  motivated  by  a  number  of  factors.  One  direct  goal  is  to  test  the 
conclusions  of  Li  and  Thurber  (1991)  regarding  the  utility  of  PmP  for  improving 
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near-regional  event  location  with  a  sparse  network,  particularly  in  terms  of 
constraining  source  depth.  The  need  for  improving  regional  event  location 
capability  is  likely  to  increase  as  concerns  regarding  nuclear  weapons  prolifera¬ 
tion  grow.  A  related  goal  is  to  reevaluate  the  events  in  Kazakhstan  studied  by 
Thurber  et  al.  (1989)  to  confirm  the  presumption  that  they  are  surface  explo¬ 
sions.  A  third  goal  is  to  provide  a  source-region  velocity  model  for  future  studies 
of  far-regional  and  teleseismic  waves  from  Kazakhstan  nuclear  explosions. 

We  study  the  regional  phases  Pg,  Pn,  and  PmP  from  a  set  of  13  events 
recorded  by  the  Natural  Resources  Defense  Council/Soviet  Academy  of  Sciences 
(NRDC/SAS)  network  in  the  Kazakhstan  region  that  operated  in  the  years 
1987  to  1988  (Fig.  1).  The  network  stations,  consisting  of  3  sets  of  three- 
component  seismometers,  produced  about  two  dozen  high-quality  multi-station 
recordings  of  explosions  in  the  epicentral  distance  range  25  to  450  km  (Thurber 
et  al.,  1989).  Our  analysis  focuses  on  data  from  stations  Bayanaul  (BAY)  and 
Karkaralinsk  (KKL),  as  station  Karasu  (KSU)  suffered  from  a  serious  site  effect 
(Berger  et  al.,  1988)  and  also  was  in  a  region  of  different  velocity  structure 
(Leith,  1987).  A  number  of  secondary  P  phases  have  been  observed  on  these 
seismograms  that  are  identified  below  as  either  PmP  or  Pg. 
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Fig.  1.  Map  of  the  Kazakhstan  region  showing  the  locations  of  stations  BAY  and  KKL  of  the 
NRDC/SAS  network  (triangles),  1987  chemical  explosions  (squares),  and  the  other  events  analyzed 
(circles). 
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Previous  Research  on  Crustal  Structure  in  Kazakhstan 

An  excellent  summary  of  studies  of  crust  and  upper  mantle  structure  in 
Kazakhstan  is  provided  by  Ryaboy  (1989).  This  work  has  shown  that  the  crust 
in  the  region  has  a  structure  with  moderate  lateral  variations  and  a  crust- 
mantle  boundary  depth  ranging  between  45  and  55  km.  Leith  (1987)  reported  a 
crustal  structure  obtained  from  a  Deep  Seismic  Sounding  (DSS)  profile  near  the 
stations  KKL  and  BAY  with  a  50-km-thick  crust  and  6  crustal  layers,  with  P 
velocities  increasing  monotonically  with  depth.  Priestley  et  al.  (1988)  used  an 
inversion  of  teleseismic  P-to-S  converted  waves  at  stations  BAY  and  KKL  to 
obtain  P  velocity  models  (assuming  a  constant  Vp/Vs  value)  generally  similar 
to  that  reported  by  Leith  (1987),  except  for  a  slightly  decreased  crustal  thick¬ 
ness  at  station  BAY,  a  broad  Moho  transition  beneath  KKL,  and  a  slight 
velocity  reversal  at  shallow  depths  beneath  both  stations  (Fig.  2).  Overall,  there 
is  evidence  for  significant  lateral  crust  and  upper  mantle  heterogeneity  in  this 
region,  though  relatively  modest  heterogeneity  in  the  immediate  vicinity  of  the 
Kazakhstan  Test  Site.  Therefore,  any  one-dimensional  crustal  model  is  only  an 
approximation  to  the  actual  structure  in  this  region.  With  this  limitation  in 
mind,  we  endeavor  to  match  reflectivity  synthetics  competed  for  a  simple 
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Fig.  2.  P- wave  velocity  models  from  Leith  (1987)  from  a  DSS  profile  in  the  Kazakhstan  region 
(solid)  and  Priestley  et  al.  (1988)  from  receiver  function  inversion  for  stations  BAY  (dotted)  and  KKL 
(two  models,  dashed  and  dot-dashed). 
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one-dimensional  crustal  model  to  observed  seismic  travel  times  and  waveforms. 
Our  main  focus  is  on  first  and  secondary  P- wave  arrivals  and  their  potential  for 
regional  event  location. 

Examination  of  Near-Regional  Data 

The  data  analyzed  in  this  study  were  recorded  on  a  set  of  surface  mounted 
Teledyne  Geotech  GS-13  short-period  seismometers,  a  set  of  Kinemetrics  SV-1 
and  SH-1  surface-mounted  intermediate-period  'eismometers,  and  a  set  of 
Geotech  54100  borehole-mounted  seismometers  installed  at  depths  of  66  and 
101  m  at  BAY  and  KKL,  respectively  (Berger  et  al.,  1988).  The  velocity  response 
of  each  of  these  instruments  (Berger  et  at.,  1988)  is  flat  from  at  least  1  to  50  Hz. 
The  noise  spectrum  for  these  stations  is  also  summarized  by  Berger  et  al. 
(1988). 

The  events  used  in  our  analysis  have  been  located  by  Thurber  et  al.  (1989) 
using  body-wave  arrival  times  and  P- wave  arrival  azimuths.  Absolute  locations 
are  known  for  the  1987  chemical  explosions  carried  out  as  part  of  the  NRDC/SAS 
verification  project  (Given  et  al.,  1990),  other  events  are  clearly  associated  with 
mines  or  quarries  observed  in  satellite  images  (Thurber  et  al.,  1989),  and  still 
others  can  be  associated  with  mines  indicated  on  regional  maps.  None  of  the 
events  used  are  thought  to  be  earthquakes.  There  is^  some  uncertainty  is  the 
epicentral  distances  of  the  events  analyzed,  so  we  place  special  emphasis  on 
those  events  whose  locations  are  most  reliable. 

An  examination  of  the  unprocessed  waveforms  at  the  NRDC/SAS  stations 
BAY  and  KKL  indicates  the  first  P  arrival  is  usually  followed  by  one  or  more 
distinct  secondary  arrivals  in  the  succeeding  0.5  to  2.5  sec  in  the  epicentral 
distance  range  between  100  and  300  km  (Fig.  3).  In  the  distance  range  under 
100  km  from  the  source,  there  are  no  clearly  identifiable  secondary  P  arrivals; 
only  the  phase  Pg  can  be  identified  (Fig.  3a).  At  about  120-km  distance  from  the 
source,  however,  the  secondary  arrival  PmP  appears  about  2.5  sec  after  the 
first  arrival.  The  amplitude  of  this  phase  relative  to  the  first  arrival  increases  in 
the  distance  range  from  120  to  170  km  from  the  source,  at  which  point  PmP 
arrives  about  2  sec  after  the  first  arrival  and  has  about  twice  the  amplitude 
(Fig.  3b).  Between  170  km  and  the  crossover  distance  of  about  240  km,  the 
phases  Pg  and  PmP  continue  converging  and  the  relative  amplitude  of  PmP 
generally  decreases,  until  by  300  km  distance  the  two  phases  are  weak  and 
indistinguishable,  and  we  have  difficulty  picking  secondary  phases  with  any 
reliability  (Fig.  3c). 

To  gain  confidence  in  our  identification  of  these  arrivals,  we  carried  out  a 
polarization  analysis  on  the  secondary  phases  to  determine  whether  they  had 
the  correct  polarization  expected  for  P  waves  from  the  event  azimuth.  The 
polarization  analysis  was  done  using  the  covariance  of  the  horizontal- 
component  seismograms,  as  described  by  Thurber  et  al.  (1989).  For  a  polarized 
signal  in  the  presence  of  noise,  the  eigenvector  corresponding  to  the  largest 
eigenvalue  of  the  covariance  matrix  for  the  signal  components  gives  the  direc¬ 
tion  of  polarization,  and  the  ratio  of  eigenvalues  measures  the  rectilinearity  of 
particle  motion  (Kanasewich,  1981).  Since  we  are  interested  just  in  the  arrival 
azimuth,  the  horizontal  component  seismograms  from  a  station  are  windowed 
and  demeaned,  and  the  2-by-2  signal  covariance  matrix  is  computed.  An  exam¬ 
ple  is  shown  in  Figure  4.  We  find  that  a  large  fraction  of  the  presumed 
secondary  arrivals  have  the  same  azimuth  as  the  first  arrival  and  show  signifi- 
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Fk;.  3.  (a  to  c)  Examples  of  observed  seismograms  at  a  range  of  epicentral  distances  A  window  of 
10  sec  is  shown  beginning  about  2  sec  before  the  first  P  arrival.  The  epicentral  distance  is  indicated 
in  the  upper  left  of  each  panel.  The  phases  Pg,  Pn,  and  PmP  are  inoicated 
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cant  particle  motion  rectilinearity.  When  comparing  the  data  to  the  synthetics, 
we  have  considered  only  the  phases  that  the  polarization  analysis  indicates  are 
from  the  same  azimuth  as  the  first  arrival.  The  polarization  analysis  and 
computed  arrival  times  and  amplitudes  of  the  synthetic  arrivals  give  us  confi¬ 
dence  in  our  ability  to  pick  the  secondary  phases  Pg  and  PmP.  Consequently, 
we  have  confined  most  of  our  attention  to  making  reliable  picks  on  these  phases, 
along  with  Pn  when  it  is  the  first  arrival. 

Crustal  Model  and  Synthetics 

After  making  initial  picks  of  the  phase  arrival  times,  we  searched  for  a  crustal 
model  that  matched  the  observed  phase  arrival  times  and  amplitude  ratios.  We 
computed  synthetic  seismograms  using  the  code  of  Mallick  and  Frazer  (1987). 
This  code  was  written  for  computing  high-frequency  (up  to  20  Hz)  synthetic 
seismograms  over  regional  distances  (0  to  1000  km)  and  a  large  number  of 
layers.  It  has  been  used  successfully  for  a  number  of  different  applications, 
including  modeling  reflection  seismograms  in  the  range  0  to  100  km  and 
modeling  high-frequency  regional  arrivals  at  distances  up  to  1400  km  (Mallick 
and  Frazer,  1987).  For  the  various  models  we  examine,  we  compute  synthetic 
seismograms  and  low-pass  filter  them  using  a  4-Hz  comer  frequency.  For 
comparison  with  the  synthetic  seismograms,  the  observed  seismograms  were 
rotated  into  the  azimuth  of  the  source,  low-pass  filtered  with  a  comer  frequency 
of  4  Hz,  and  then  plotted  for  10  sec  around  the  first  P- wave  arrival. 

To  improve  our  crustal  model,  we  first  computed  arrival  times  of  Pn,  PmP, 
and  Pg  to  compare  with  the  data.  The  crustal  models  were  varied  systemati- 
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Fig.  4.  Polarization  analysis  for  5  sec  of  an  example  event  showing  estimated  arrival  azimuth 
(thin  solid  line)  and  calculated  particle  motion  rectilinearity  (broad  shaded  line)  for  the  Prt,  Pg.  and 
PmP  phases.  Arrival  azimuth  ranges  are  plotted  for  the  range  ±  90",  while  rectilinearity  is  plotted 
in  the  upper  portion  with  an  arbitrary  scale.  This  event  is  the  254  km  distance  record  in  Figure  3 
(chemex  1,  872450700,  at  KKL). 


cally  to  see  the  effects  on  the  synthetic  seismograms.  We  started  with  the  Leith 
(1987)  crustal  model  and  tested  variations  in  the  upper-  and  mid-crustal  veloc¬ 
ity  structure  until  we  obtained  good  agreement  with  the  Pg  first  arrivals  in  the 
distance  range  less  than  100  km  from  the  events.  Pg  and  PmP  arrival  times  in 
the  distance  range  100  to  245  km  helped  determine  the  lower  crustal  velocities. 
Near  the  crossover  distance  of  about  240  km,  we  find  considerable  waveform 
complexity  in  both  the  data  and  synthetics.  We  modeled  the  Pn  first  arrivals  in 
the  distance  range  between  240  and  300  km  to  ascertain  the  transition  velocity 
between  40  and  50-km  depth.  A  transition  velocity  of  6.95  to  7.05  km/sec  was 
needed  in  this  depth  range  to  model  the  arrival  times  of  the  data.  The  data 
require  a  sharp  transition  at  the  crust-mantle  boundary  at  50-km  depth  to 
match  the  observed  difference  in  the  travel  times  of  the  different  phases.  We 
constrained  the  upper-mantle  velocity  in  the  depth  range  between  50  and  85  km 
using  the  Pn  arrival  times  in  the  distance  range  between  250  and  350  km.  A 
velocity  of  8.25  km/sec  in  the  upper  mantle  between  55-  and  85-km  depth  gave 
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the  best  fit  to  the  arrival  time  of  Pn  and  the  observed  difference  in  the  arrival 
times  of  Pn  and  PmP  in  this  distance  range.  Beyond  about  260  km,  the  largest 
amplitude  secondary  arrival  has  the  arrival  time  expected  of  PmP\  however, 
the  synthetic  PmP  considerably  exceeds  the  observed  amplitude  of  PmP  in  the 
data. 

Our  final  velocity  model  is  presented  in  Figure  5.  It  has  an  upper-crustal 
velocity  increasing  fairly  rapidly  from  4.5  km/sec  near  the  surface  to  6.5 
km/sec  at  15-km  depth,  then  increasing  more  slowly  to  7.05  km/sec  at  50-km 
depth.  The  crust-mantle  transition  zone  ranges  from  8.05  km/sec  at  50-km 
depth  to  8.25  km/sec  at  55-km  depth.  The  upper  mantle  has  a  P  velocity  of 
about  8.25  km /sec  between  55-  and  85-km  depth.  Below  this  depth,  the  data  do 
not  constrain  the  velocity  structure,  and  we  have  used  the  structure  obtained  by 
Goldstein  et  al.  (1992).  We  have  assumed  a  constant  value  for  Vp/Vs  of  1.73 
throughout  the  model.  Overall,  our  final  model  is  similar  to  that  obtained  from 
the  deep  seismic  soundings  (dashed  lines  in  Fig.  5). 

The  fit  of  the  data  to  the  synthetics  is  quite  good  in  the  distance  range 
between  100  and  280  km  (Fig.  6)  and  is  superior  to  the  Leith  model  (Fig.  7). 
This  model  matches  the  Pg,  PmP,  and  Pn  arrival  times  within  about  0.3  sec  for 
most  of  the  phases  seen  on  the  seismograms.  We  get  very  good  agreement  with 
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Fin.  5.  Final  velocity  structure  model  for  the  Kazakhstan  region  (solid  line)  compared  to  the 
initial  model  from  Lcitn  (1987)  (dashed  line). 
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Fig.  6.  (a  and  b)  Record  sections  of  (a)  synthetic  seismograms  at  100,  150,  200,  250,  and  300  km 
distances  and  (b)  observed  data  at  comparable  distances  (871460833,  872450927,  871351035, 
872450700,  and  872440952),  using  a  reducing  velocity  of  7.5  km/sec. 
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the  Pg  first  arrival  times  in  the  distance  range  less  than  100  km;  however,  we 
find  that  the  initial  coda  of  the  data  are  larger  in  amplitude  than  the  synthetics. 
We  can  attribute  this  to  wave  scattering  and  conversion  not  modeled  by  the 
reflectivity  synthetics.  Between  100  and  200  km,  the  first  arrival  times  and  the 
overall  waveform  envelopes  show  good  agreement;  however,  the  arrival  time 
difference  between  Pg  and  PmP  is  greater  in  the  synthetics  than  in  the  data  by 
between  0.2  and  0.7  sec  for  the  distance  range  under  150  km,  and  less  in  the 
synthetics  than  in  the  data  by  between  0.3  and  0.9  sec  for  the  distance  range 
beyond  150  km.  The  model  shows  the  best  agreement  in  the  distance  range 
between  180  and  240  km,  where  the  synthetics  match  both  the  arrival  times 
and  amplitudes  of  Pg  and  PmP  within  0.3  sec  and  50%,  respectively.  Near  the 
crossover  distance  of  240  km,  we  find  that  the  Pg,  Pn,  and  PmP  phases  all 
arrive  in  the  first  1.5  sec  of  the  record  on  both  the  synthetic  and  the  data 
records.  Between  250  and  280  km,  the  synthetic  matches  the  arrival  time  of  Pn 
and  PmP  within  0.3  sec;  however,  the  Pn  amplitude  of  the  synthetic  underesti¬ 
mates  that  of  the  data  while  the  synthetic  PmP  amplitude  exceeds  that  of  the 
data,  each  by  about  a  factor  of  2.  Beyond  300  km  we  have  difficulty  matching 
the  synthetics  to  the  data:  every  realistic  model  that  matches  the  Pn  and  PmP 
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Fig.  7.  (a  to  e)  Comparison  of  observed  seismograms  (top)  with  reflectivity  synthetic  seismograms 
at  105,  156,  197,  257,  and  301  km  distance  for  our  final  model  (middle)  and  the  Leith  model 
(bottom).  The  amplitude  scale  is  arbitrary. 
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Fig.  7.  ( Continued ). 

arrival  times  in  the  distance  range  under  300  km  gives  a  PmP  amplitude  that 
is  a  factor  of  2  greater  than  that  observed  in  the  data  beyond  300  km.  We 
believe  that  lateral  heterogeneity  probably  renders  the  assumption  of  a  simple 
layered  model  questionable  beyond  300  km.  Therefore,  we  believe  that  using  a 
single-layered  crustal  model,  as  required  for  the  reflectivity  synthetics,  is  proba¬ 
bly  inappropriate  for  events  in  the  distance  range  beyond  300  km. 

Event  Relocation  Using  the  Phases  Pg,  Pn,  and  PmP 

Previous  empirical  and  theoretical  work  on  regional  event  location  in 
Kazakhstan  (Thurber  et  al.,  1989;  Li  and  Thurber,  1991)  has  investigated  the 
usefulness  of  secondary  P  arrivals  for  regional  event  location.  With  an  improved 
crustal  model  for  Kazakhstan  (Fig.  5)  and  a  determination  of  Pg  and  PmP 
secondary  arrivals  in  the  data  supported  by  synthetic  seismogram  modeling 
(Figs.  6  and  7),  we  relocate  the  regional  events  in  Kazakhstan.  Event  locations 
were  computed  in  the  same  manner  as  described  by  Thurber  et  al.  (1989), 
incorporating  arrival  azimuth  information  and  utilizing  the  location  algorithm 
TTAZLOC  (Bratt  and  Bache,  1988). 

To  evaluate  the  utility  of  PmP  for  event  location,  13  of  the  events  shown  in 
Figure  1  for  which  PmP  could  be  identified  were  relocated  using  the  new 
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Fig.  7.  ( Continued ). 


velocity  model,  both  with  and  without  the  phase  PmP.  We  identified  PmP  on 
nearly  half  the  seismograms  (see  Fig.  3  for  examples),  with  estimated  reading 
uncertainties  of  0.05  to  0.25  sec.  Initally,  the  events  were  relocated  with  focal 
depth  fixed  at  0  km,  as  had  been  done  by  Thurber  et  al.  (1989)  for  all  but  the 
three  chemical  explosions  (Given  et  al.,  1990).  The  event  locations  showed  little 
change,  with  differences  averaging  less  than  4  km  in  latitude  and  5  km  in 
longitude.  The  mean  RMS  residual  for  all  PmP  observations  was  about  0.35  sec, 
compared  to  0.25  and  0.50  sec  for  all  Pn  and  Sn  observations,  respectively. 
Thus,  the  PmP  arrival  times  were  of  a  quality  comparable  to  the  previous  Pn 
and  Sn  data.  The  estimated  location  uncertainties  decreased  about  10%  when 
PmP  observations  were  added,  again  indicating  that  the  quality  of  the  PmP 
data  is  comparable  to  that  of  the  other  phases.  Most  of  this  decrease  can  be 
attributed  to  the  increase  number  of  degrees  of  freedom  due  to  adding  PmP 
without  an  increase  in  data  misfit. 

The  theoretical  study  of  Li  and  Thurber  (1991)  indicates  that  PmP  arrival 
times  can  provide  significant  constraint  on  source  depth  for  regional  events 
recorded  by  a  sparse  network.  Therefore,  we  recomputed  the  event  locations 
with  focal  depth  unconstrained  both  with  and  without  the  PmP  observations. 
The  starting  value  for  depth  was  0  km,  and  negative  depth  values  (i.e., 
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“airquakes”)  were  prohibited  by  the  algorithm.  The  RMS  residuals  from  the 
constrained  focal  depth  location  results  above  were  adopted  as  a  priori  values 
for  the  data  variances.  As  in  previous  studies  (Jordan  and  Sverdrup,  1981;  Bratt 
and  Bache,  1988;  Thurber  et  al.,  1989),  a  priori  information  was  given  a 
K-weight  of  8. 

The  location  results  with  focal  depth  free  are  indicated  in  Table  1.  The 
epicenter  locations  were  again  similar  with  or  without  PmP;  all  events  had 
differences  in  epicenter  less  than  5  km.  Epicentral  uncertainties  were  generally 
reduced  with  PmP,  but  not  substantially.  The  situation  for  focal  depth  was 
quite  different,  however.  Although  computed  depth  differences  were  not  signifi¬ 
cant,  nearly  half  of  the  events  had  their  depth  uncertainty  estimate  reduced  by 
a  factor  of  1.5  or  greater  when  PmP  observations  were  included.  Thus,  the  PmP 
observations  do  provide  useful  constraint  on  source  depths  at  these  near- 
regional  distances. 

An  important  issue  implicit  in  this  discussion  is  the  confidence  that  these 
events  are  explosions  and  not  earthquakes.  In  Table  1,  independent  information 
on  or  associations  of  these  explosions  with  mines  or  quarries  is  indicated.  Two  of 
the  events  are  the  1987  chemical  explosions,  and  suspected  source  areas  for 
seven  other  events  had  been  identified  previously  by  Thurber  et  al.  (1989)  from 
satellite  images  (mines  at  Ekibastuz  and  Balkash).  Of  these,  events  871351035, 
871410916,  and  871460833  were  also  shown  to  have  spectral  modulation  typical 
of  ripple-fired  blasts  by  Hedlin  et  al.  (1989),  In  addition,  one  event  (871430849) 
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is  thought  to  have  been  a  blast  at  a  known  quarry  in  the  town  of  Karagayly  (H. 
Given,  personal  comm.).  Of  the  remaining  four  events,  event  871351035  is 
located  within  4  km  of  a  mapped  mine  at  the  town  of  Yuzhnyy  on  the  1983  ONC 
navigational  chart  for  the  region,  while  event  871410916  is  located  within 
10  km  of  two  mapped  mines  at  the  town  of  Molodezhnoye.  Of  the  other  two 
events,  871340936  occurred  in  the  Karaganda  area,  known  for  mining,  at  a  time 
of  day  typical  for  mine  blasts.  Only  event  871460531  is  “suspicious”  due  to  its 
time  of  occurrence,  although  the  location  results  are  not  inconsistent  with  a 
surface  focus. 

Conclusions 

Using  a  reflectivity  synthetic  seismogram  code,  we  have  modeled  primary  and 
secondary  P  phases  for  a  data  set  from  near-regional  events  recorded  in  1987 
and  1988  by  the  stations  BAY  and  KKL  of  the  former  NRDC/SAS  network.  An 
analysis  of  wave  polarization  was  used  to  help  identify  the  secondary  phases, 
primarily  PmP.  A  new  layered  crustal  model  for  the  region  was  developed  to 
improve  the  fit  to  the  arrival  times  and  waveforms  of  these  phases.  We  can 
match  a  Pg  first  arrival  and  a  PmP  second  arrival  on  most  of  the  seismograms 
in  the  distance  range  between  100  km  and  the  crossover  distance  of  240  km 
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with  an  arrival  time  misfit  on  the  order  of  0.5  sec,  and  the  relative  amplitudes 
are  also  matched  within  50%.  Beyond  the  crossover  distance,  we  are  able  to 
model  Pn  and  PmP  arrival  times  with  comparable  fit,  but  we  have  difficulty  in 
matching  the  observed  amplitude  of  the  phase  PmP  beyond  300  km,  presum¬ 
ably  due  to  effects  of  lateral  heterogeneity.  A  secondary  Pg  phase  is  only 
observed  in  the  range  just  beyond  the  crossover  distance. 

Our  crustal  model  for  the  Kazakhstan  region  has  a  4.5  to  6.5  km/sec  upper 
crust  (to  15-km  depth),  a  middle  and  lower  crust  that  increases  in  velocity  from 
6.5  to  7.05  km/sec,  and  a  crust-mantle  boundary  at  about  50-km  depth. 
Compared  to  the  model  of  Leith  (1987)  derived  from  DSS  studies,  our  upper 
crust  has  a  slightly  lower  velocity  and  our  lower  crust  has  a  slightly  higher 
velocity,  similar  to  the  results  of  Priestley  et  al.  (1988).  A  relatively  sharp 
crust-mantle  boundary  (about  5-km  thick)  was  needed  to  match  the  observed 
difference  in  the  Pn  and  PmP  arrival  times.  This  is  in  contrast  to  the  receiver 
function  results  of  Priestley  et  al.  (1988)  for  the  structure  beneath  station  KKL, 
which  suggested  a  10-  to  15-km-thick  transitional  Moho.  The  uppermost  mantle 
has  a  velocity  of  8.05  km/sec  between  50-  and  55-km  depth,  increasing  to  8.25 
km/sec  between  55-  and  85-km  depth.  The  near-regional  data  do  not  constrain 
the  structure  at  greater  depths. 

The  new  velocity  structure  and  PmP  observation^  obtained  from  this  study 
were  used  to  relocate  events  in  the  distance  range  up  to  about  300  km  from  the 
stations.  It  was  found  that  adding  the  phase  PmP  reduces  the  uncertainty  for 
event  depth  in  all  cases,  with  relatively  little  effect  on  epicenter  uncertainty. 
None  of  the  events  can  be  conclusively  demonstrated  to  be  earthquakes,  based 
on  their  focal  depths  and  associations  with  active  mines  and  quarries,  though 
one  event  cannot  be  confidently  classified  as  an  explosion. 
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2.4.2.  Modified  generalized  ray  modeling  of  regional 
seismograms  in  Central  Asia 

Introduction 

The  demise  of  the  Soviet  Union  coupled  with  increased  concerns  of  nuclear  weapons 
proliferation  monitoring  is  resulting  in  more  research  effort  in  the  area  of  seismic  wave  propagation 
at  regional  distances  (500  to  2000  km).  This  distance  range  is  difficult  to  model  accurately  because 
it  is  too  short  to  permit  a  teleseismic-like  approach  in  which  only  the  rece.ver-site  influence  on  the 
waveform  is  of  concern,  and  it  is  too  long  for  a  single  source-region  velocity  model  to  be  valid. 
Thus,  waveform  modeling  methods  that  incorporate  lateral  heterogeneity  are  required. 

A  few  researchers  have  attempted  to  model  waveforms  and  crustal  structure  in  the  C.I.S. 
and  Kazakhstan  regions  using  single  structure  reflectivity  or  generalized  ray  models,  with  varying 
levels  of  success.  Goldstein  et  al.  (1992)  modeled  a  number  of  nuclear  explosions  at  station  ARU 
using  a  reflectivity  code.  Their  research  gave  a  model  of  the  mantle  between  50  and  400  km  depth 
in  the  region  between  the  Kazakhstan  test  site  and  station  ARU.  This  model  showed  a  low  velocity 
zone  in  the  depth  range  between  1 50  and  200  km.  Burdick  et  al.  ( 1992)  modeled  data  recorded  at 
about  10  stations  in  Kazakhstan  and  the  C.I.S.  at  distances  between  1500  and  about  5000  km  for 
events  recorded  in  the  years  between  1971  and  1989,  using  pre-existing  velocity  structures  and  the 
generalized  ray  method.  They  made  no  attempt  to  develop  new,  improved  models.  Both  of  these 
studies  share  two  similar  problems;  first  they  require  a  common  model  at  both  the  source  and 
receiver,  which  does  not  allow  for  the  effect  of  lateral  heterogeneity  in  this  region,  and  second  an 
inability  to  modei  a  large  amount  of  the  observed  complexity  of  the  data,  especially  the  secondary 
arrivals  in  the  time  frame  of  5  to  20  seconds  behind  the  first  arrival. 

In  a  study  of  the  receiver  structure  under  the  sites  ARU,  GAR,  KIV,  and  OBN,  Riviere- 
Barbier  et  aJ.  (1992)  used  teieseisrric  P  waveforms  to  model  the  velocity  structure  under  these 
stations.  They  produced  models  which  match  the  teleseismic  data  fairly  well,  deriving  significantly 
different  receiver  structures  for  each  station.  However,  their  results  do  not  match  those  of  Burdick 
et  al.  (1992)  or  Goldstein  et  al.  ( 1992).  At  station  WMQ,  M'  .igino  and  Ebel  (1992)  modeled 
teleseismic  data  to  obtain  a  preliminary  receiver  crustal  structure;  however,  they  do  not  resolve  the 
details  of  the  upper  mantle  structure  below  90  km  depth. 

In  this  section,  we  present  crustal  and  upper  mantle  models  under  the  stations  BRV,  FRU, 
NVS,  and  WMQ  (Figure  1)  which  for  FRU  and  NVS  are  the  first  known  models  at  this  site,  and 
which  for  BRV  and  WMQ  are  the  first  models  of  the  upper  mantle.  They  are  obtained  from  a 
modified  generalized  ray  synthetic  seismogram  computation  using  a  split  source-receiver  model 
which  computes  approximate  synthetic  seismograms  for  the  first  10  or  20  seconds.  Our  results 
give  a  better  fit  to  the  data  by  nearly  a  second  in  first  and  second  arrival  times  over  results  obtained 
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using  a  common  source-receiver  model.  The  source  region  structure  in  the  depth  range  between  0 
and  90  km  is  the  Kazakhstan  model  obtained  from  reflectivity  synthetic  seismograms  of  chemical 
explosions  determined  by  Quin  and  Thurber  (1992).  The  mantle  structure  in  the  depth  range  from 
90  to  200  km  was  obtained  by  matching  the  regional  synthetics  in  the  distance  range  of  600  to 
1000  km.  The  upper  mantle  models  show  evidence  for  a  moderate  amount  of  lateral  heterogeneity 
in  this  region;  most,  however,  show  an  increase  from  8.25  km/sec  at  about  50  km  depth  through 
8.37  km/sec  at  70  km  depth  to  about  8.55  km/sec  at  100  km  depth. 

Generalized  Ray  Model 

The  modeling  procedure  used  in  this  section  uses  the  generalized  ray  method  of 
Helmberger  (1972)  but  with  a  significant  modification.  The  original  generalized  ray  method 
assumes  that  the  source  and  receiver  have  the  same  velocity  model.  However,  regional  refraction 
data  indicates  that  the  C.I.S.  and  Kazakhstan  region  exhibit  considerable  crust  and  mantle  lateral 
heterogeneity  (Ryaboy,  1989).  An  examination  of  the  Cagniard  integration  method  indicates  that 
the  amplitude  and  phase  of  the  arriving  waves  can  be  broken  down  into  separate  functions  of  the 
source  and  receiver  models,  with  the  down-going  wave  propagating  through  the  source  region  and 
the  up-going  wave  propagating  through  the  receiver  region  for  the  phases  Pn,  PmP  and  reflected 
waves  which  propagate  through  the  mantle  (Figure  2).  The  Laplace  transform,  s,  of  the 
generalized  ray  response  is  given  by  (Aki  and  Richards,  1980) 

(1)  P(r,z,s)=2s/JtRsP0(s)Imj^  Ko(spr)[PR(p)]exp(-s[SUM(p)])dp 
where 

(2)  PR(p)=(PP)ls(PP)2r  •  (PP)(n-l)s  (PP)n,(PP)(n-Dr-  •  (PP)2r(PP)lr 

where  PPng  and  PPnr  indicate  the  reflection  and  refraction  coefficients  of  down-going  (source)  and 
up-going  (receiver)  waves,  respectively,  and  where 

(3)  SUM(p)  =  (Thi,  -  ds)^js  +  Th2S  •  •  +  Th(n.i)S  £<n  i  >«,  +  Thns  +  Th(n.1)r  ^(„.i)r 

•  •  •  Th2r  ^2r  +(Thir  -dr)^ir 

where  Thj  is  the  thickness  of  each  layer  and  d5  and  dr  are  the  som .  *■  and  receiver  depths, 

4i  =  (aj2  -  p 2),/2  is  the  ray  parameter  weighting  factor  for  each  layer,  Ro  is  the  distance  from  the 
source  to  the  receiver,  Pq(s)  is  the  Laplace  transform  of  the  initial  pressure  pulse,  and  Ko(spr)  is  a 
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modified  Bessel  function.  The  Cagniard  path  fen-  inversion  of  equation  1  is  now  the  solution  of 
p  =  p  (t)  where  x  =  pr  +  sum(p). 

The  split  source  receiver  method  entails  solving  for  the  functions  PR(p)  and  SUM(p). 
These  equations  can  be  split  into  two  parts:  a  downward  part  PPisPP2SPP3s—  and  an  upward  part 
. . .  .PP3r  PP2r  PPir  (note  that  for  the  case  of  converted  waves  we  would  replace  the  terms  PP  by 
PS,  SS  or  SP).  For  SUM(p)  we  need  to  solve  the  equation  p  =  p(x)  for  t  numerically  in  order  to 
find  the  Cagniard  integration  path  for  each  ray.  This  solution  begins  with  an  initial  ray  parameter 
on  the  real  axis  obtained  from  the  P  wave  slowness  of  the  deepest  layer  in  which  that  ray 
penetrates.  Subsequent  values  of  the  complex  ray  parameter  p  for  each  time  point  on  the 
seismogram  are  obtained  from  an  iterative  method  in  which  the  initial  guess  of  the  complex  ray 
parameter  arriving  at  that  time  step  is  obtained  from  the  last  ray  parameter  and  the  direction  of  the 
ray  parameter  contour  at  the  last  time  step. 

In  terms  of  modification  of  the  generalized  ray  program,  the  subroutine  which  determines 
the  Cagniard  travel  time  is  modified  to  incorporate  the  split  source  receiver  structure,  with  a  slight 
change  in  the  ray  parameter  computation  routine  to  account  for  rays  which  do  not  converge.  In 
practice  the  solution  to  this  equation  can  be  found  for  about  99%  of  all  possible  rays  arriving  in  the 
first  20  seconds  after  the  first  P  wave,  including  all  of  the  most  important  rays;  the  only  rays  for 
which  a  ray  parameter  path  cannot  be  found  are  for  a  few  late  arriving  highly  converted  phases. 
Therefore,  we  believe  our  method  is  an  effective  approximation  procedure  for  the  split  source- 
receiver  case. 

Our  model  treats  the  source-receiver  boundary  as  being  exactly  at  the  mid-point  of  the 
source-receiver  distance.  In  reality,  this  boundary  could  be  located  anywhere;  and,  in  theory,  the 
wave  integral  could  be  split  up  so  that  multiply  reflected  phases  could  be  treated  as  if  they 
propagated  through  more  than  two  distinct  regions  with  different  phase  factors  and  different 
refraction  coefficients  for  each  set  of  layers,  as  long  as  a  solution  for  the  ray  parameter  equation 
can  be  found.  However,  we  know  from  the  work  of  Quin  and  Thurber  (1992;  section  2.4.1)  that 
the  Kazakhstan  region  does  not  show  significant  heterogeneity  in  the  distance  range  less  than  about 
300  km  from  the  source,  and  that  most  lateral  heterogeneity  occurs  in  the  crustal  region,  where  the 
model  of  waves  traveling  downwards  from  the  source  and  upwards  to  the  receiver  is  most 
appropriate.  Also,  the  upper  mantle  in  this  region  shows  only  limited  heterogeneity;  computing  the 
mantle  velocity  for  head  waves  using  the  average  mantle  velocity  gives  an  incorrect  velocity  by  at 
most  0.03  km/sec,  well  within  the  accuracy  of  the  known  arrival  times.  In  addition,  surface 
reflection  phases  2PmP  and  2Pn  can  be  considered  to  propagate  through  a  source  region,  bounce 
off  the  free  surface,  and  propagate  through  a  receiver  region  of  different  crust  and  upper  mantle 
structure.  Therefore,  we  can  use  the  same  split  source-receiver  approximation  for  these  phases  as 
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well.  This  means  that  we  can  model  almost  all  phases  arriving  in  the  first  20  seconds  in  the 
distance  ranges  of  interest  for  these  seismograms  using  the  split  source-receiver  model. 

The  synthetic  seismograms  are  sampled  at  20  Hz  and  inverse  Fourier  transformed  at  a 
sampling  length  of  1024  points,  and  then  low  pass  filtered  with  a  comer  frequency  of  2  Hz.  This 
enables  us  to  get  a  good  sampling  of  the  important  phases,  and  to  filter  out  high  frequency  noise. 
We  incorporate  attenuation  using  the  standard  constant  Q  attenuation  model  with  a  Q  value  of  500, 
a  value  typical  of  crustal  P  waves  in  this  region  (Sereno,  1990).  This  gives  us  a  reasonable  value 
for  attenuation  in  this  region,  eliminating  unphysical  high  frequency  arrivals.  Using  this  method 
we  can  reproduce  most  of  the  important  first  arriving  phases  seen  in  the  time  frame  of  the  first  20 
seconds  of  data. 

To  test  the  accuracy  of  the  program  we  compared  synthetic  seismograms  from  a  single 
model  with  those  computed  by  the  separate  source-receiver  model  method  for  the  case  of  identical 
crustal  structures  and  find  that  the  two  methods  give  identical  results.  For  the  case  where  the  two 
models  are  different  we  have  found  that  the  split  model  enables  us  to  produce  a  much  better  fit  to 
the  arrival  times  and  phases  than  the  single  structure  model  (Figure  3).  This  figure  shows  both 
BRV  and  WMQ  waveforms  computed  at  a  distance  of  670  and  800  km.  Our  split  model  gives  a 
better  fit  to  the  data  than  the  single  model.  In  particular,  we  are  able  to  accurately  match  both  the 
first  and  second  arrival  times  and  the  overall  waveform  envelopes  better  with  our  approach.  We 
believe  our  approach  more  effectively  incorporates  the  known  cmstal  and  upper  mantle  lateral 
complexity  as  determined  by  refraction  and  teleseismic  data  than  efforts  based  on  single  structure 
models. 


Regional  Data 

Data  from  a  number  of  different  sources  were  utilized  in  this  study.  Tables  1  and  2 
summarize  the  stations  and  event  data  used  in  this  study.  The  quality  and  type  of  data  used  varied 
from  station  to  station.  For  stations  FRU  and  NVS  the  data  consist  of  intermediate  period  vertical 
hand-digitized  data  from  explosions  in  the  Kazakhstan  and  Konystan  regions  in  1965-1989;  at 
stations  BRV  and  WMQ  the  data  consist  of  high  quality  broad  band  3-component  recordings  made 
in  1988.  All  of  the  data  at  stations  WMQ  and  BRV  are  usable;  however,  a  large  number  of  the 
waveforms  at  stations  FRU  and  NVS  are  incomplete  or  inaccurate. 

We  have  taken  the  usable  waveforms  from  sets  of  events  with  epicenters  within  about  3  km 
of  each  other  as  determined  from  satellite  images  and  teleseismic  relocations  (Lilwall  and  Farthing, 
1990,  Thurber  et  al.,  1992),  deconvolved  the  station  responses  and  stacked  them  to  obtain  3  or  4 
different  modelable  waveforms  at  each  station.  A  large  number  of  these  waveforms  are  clipped 
and  inaccurately  timed;  however,  about  a  dozen  or  so  were  usable.  These  waveforms  span  four 
distances  in  the  range  621-676  km  at  stations  NVS  and  813-847  km  at  station  FRU.  At  stations 
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BRV  and  WMQ,  we  low  pass  filtered  with  a  comer  frequency  of  2  Hz  and  rotated  the  broadband 
data  from  several  nuclear  explosions  recorded  in  the  years  1972- 1989  to  get  good  quality 
waveforms  in  the  distances  between  683-697  km  at  BRV  and  950-962  km  at  WMQ.  After 
comparing  both  stacked  and  unstacked  data  at  BRV  and  WMQ  against  the  synthetics,  we  decided 
to  model  the  unstacked  data,  due  to  the  high  quality  of  the  data  and  the  high  accuracy  with  which 
the  explosion  locations  after  1986  were  known  from  satellite  and  teleseismic  data  (Thurber  et  al., 
1992;  section  2.4.3). 


Explosion  Source  Model 

The  manner  in  which  the  characteristics  of  the  nuclear  explosion  source  are  quantified  has 
been  studied  in  considerable  detail  in  both  the  near  and  far  field.  In  general,  the  details  of  the 
manner  in  which  the  source  is  described  depends  greatly  upon  the  distance  range  at  which  it  is 
observed;  in  the  far  field,  the  explosion  source  appears  simply  an  as  isotropic  point  force  right 
below  the  surface,  while  in  the  near  field  the  source  appears  as  a  complex  nonlinear  shock  wave. 
Day  (1983)  and  Von  Seggem  (1988)  have  studied  the  effect  of  source  depth  and  spallation  between 
100  and  1000  km  from  the  explosion  source.  In  particular,  they  have  concluded  that  both  source 
depth  and  spallation  are  important  in  the  generation  of  the  regional  phase  Lg,  especially  since  this 
phase  cannot  be  successfully  modeled  using  simple  planar  models  and  linear  point  explosion 
terms. 

Their  computation  of  the  phases  Pn  and  pPn,  the  near  surface  reflected  phase,  indicate  that 
to  first  order  the  near  surface  velocity  and  the  source  depth  arc  the  most  important  parameters 
controlling  the  generation  of  these  phases  in  the  mid-regional  case.  Nonlinear  effects  appear  as 
second  order  effects  which  cannot  be  distinguished  in  the  data  from  models  in  which  the  timing 
and  amplitude  of  these  phases  are  modeled  using  an  elastodynamic  representation  theorem.  In 
addition,  our  generalized  ray  method  does  not  allow  for  nonlinear  source  effects  to  be  incorporated 
in  our  models.  Therefore  in  this  report  we  have  ignored  nonlinear  near  source  effects  and  simply 
approximated  the  explosions  as  point  sources  with  an  isotropic  moment  tensor  at  an  average  depth 
of  0.5  km,  a  depth  chosen  as  the  average  obtained  from  empirical  magnitude  estimates  of  depth 
(Jih  and  Wagner,  1991),  and  included  the  surface  reflection  phases  pPg,  and  pPn.  For  the 
explosion  source  time  function  we  use  the  standard  explosion  model  of  Mueller  and  Murphy 
(1971)  which  consists  of  a  sinusoid  multiplied  by  an  exponential  term,  which  simulates  the  effect 
of  a  sharp  outward  pulse,  followed  by  an  inward  elastic  rebound. 

Modeling  Procedure 

The  modeling  procedure  consisted  of  a  two  stage  process.  In  stage  1,  we  modeled  the 
arrival  time  and  phase  of  the  first  arrivals  using  the  phase  Pn,  which  a  travel  time  analysis  indicates 
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was  the  first  arriving  phase  in  the  region  between  600  and  1000  km.  We  have  an  adequate  model 
of  the  crust  and  upper  mantle  structure  in  the  Kazakhstan  source  region  from  the  earlier  modeling 
work  of  Quin  and  Thurber  ( 1992).  Our  Kazakhstan  crust  and  upper  mantle  structure  (Figure  5) 
has  a  low  velocity  crust  of  5.0  km/sec  in  the  upper  10  km,  a  middle  crustal  layer  of  about  6.95 
km/sec  and  an  upper  mantle  velocity  of  8.25  km/sec.  Between  70  and  95  km  depth  the  source 
models  have  a  velocity  of  8.35  km.  We  have  adopted  this  crust  and  upper  mantle  model  as  the 
source  model  for  our  study;  all  of  the  rays  propagate  downward  through  this  crust  and  upper 
mantle  model  in  the  region  between  0  and  70  km  depth.  With  the  source  structure  fixed,  we  use 
the  Pn  arrival  times  and,  starting  with  crustal  models  based  on  earlier  work,  determine  the  best 
fitting  receiver  crust  and  uppermost  mantle  structures  in  each  region.  This  fixes  the  receiver  model 
down  to  about  50  km. 

With  the  crustal  structure  determined,  we  determine  the  upper  mantle  structure  by  modeling 
the  distinct  secondary  phases  observed  in  the  data  arriving  up  to  10  seconds  after  the  Pn  phase. 

Our  investigation  indicates  that  these  second  arrivals  are  generally  influenced  by  the  structure  in  the 
upper  mantle;  in  particular,  the  arrivals  in  the  first  5  seconds  after  the  Pn  arrival  are  affected  by 
structure  in  the  50  to  90  km  depth  range,  while  arrivals  in  the  5-15  sec  time  span  are  affected  by  the 
90  to  200  km  depth  range.  We  believe  our  models  resolve  details  of  the  crust  and  upper  mantle 
structure  in  this  region  down  to  a  depth  of  about  200  km  for  events  in  the  distance  range  up  to 
1000  km  from  the  source. 

We  allowed  the  source  mantle  velocity  in  the  region  between  70  and  150  km  depth  to  be  a 
variable  in  our  study.  This  at  first  may  appear  inconsistent;  however  an  examination  of  the  path  of 
the  waves  indicates  that  they  are  separated  at  this  depth  by  as  much  as  500-600  km,  which  means 
that  the  waves  propagate  through  different  tectonic  regions  towards  each  station.  We  allow  the 
velocity  structure  to  vary  by  about  0.25  km/sec  between  the  various  models  in  this  region, 
consistent  with  the  observed  lateral  heterogeneity  in  this  region  (Riviere-Barbier  et  al.,  1992);  the 
source  model  at  station  BRV  shows  a  velocity  of  8.25  km  ./sec  in  the  region  between  70  and  90 
km,  and  a  low  velocity  zone  in  the  region  between  1 10  and  1 30  km  while  the  other  models  didn't. 

For  waves  propagating  to  stations  NVS  and  WMQ  the  best  fitting  source  model  had  a 
velocity  of  8.35  km./sec  at  70  and  90  km  depth,  and  a  velocity  of  8.50  km/sec  in  the  region 
between  90  and  1 10  km  depth.  At  station  FRU,  we  find  that  the  best  fitting  source  model  had  a 
high  velocity  region  of  8.50  km/sec  in  the  region  between  90  and  1 30  km  depth.  In  the  region 
between  1 10  and  130  km,  NVS  and  FRU  have  a  velocity  of  about  8.55  km/sec,  while  all  the 
models  have  a  velocity  of  8.60  km/sec  below  1 30  km  depth  down  to  200  km  depth.  This  variation 
of  the  source  crustal  and  upper  mantle  model  at  depth  allowed  us  to  get  a  better  fit  to  the  synthetic 
seismograms. 
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About  a  dozen  trial  and  error  modeling  attempts  were  required  at  stations  FRU  and  NVS, 
for  which  we  did  not  have  good  starting  models.  At  stations  WMQ  and  BRV,  for  which  we  have 
useful  starting  models  obtained  from  teleseismic  data,  only  about  four  or  five  trials  were  needed. 
The  fit  of  the  synthetics  to  the  data  varies  from  station  to  station  (Figure  4).  At  station  NVS,  we 
match  the  first  arrivals  and  the  timing  of  the  second  and  third  arrivals  fairly  well;  however,  we  have 
difficulty  in  matching  the  amplitude  of  the  second  arrival  about  2  seconds  after  the  first  arrival.  At 
station  FRU,  the  first  arrival  of  the  synthetic  is  about  a  half  second  early,  and  we  find  that  the  third 
synthetic  arrival  about  6  seconds  after  the  first  P  wave  come  in  too  strongly,  regardless  of  what 
mantle  model  we  utilize.  We  do  match  the  timing  and  amplitude  of  the  second  arrivals  in  the  period 
about  2-3  seconds  after  the  first  arrival. 

At  stations  BRV  and  WMQ,  we  are  able  to  model  both  the  vertical  and  rotated  radial 
waveforms.  At  station  BRV,  we  are  able  to  match  both  the  first  arrival  time  of  the  phase  Pn  and 
the  overall  waveform  envelope  in  the  following  10  seconds;  however,  our  model  somewhat 
overestimates  the  first  arrival  amplitude  compared  to  the  coda.  At  station  WMQ,  we  match  the  first 
arrivals  and  the  overall  waveform  envelope.  We  do  not  get  a  good  match  to  the  second  arrivals  in 
the  next  few  seconds  after  the  first  arrival  on  the  vertical  component,  although  on  the  radial 
component  we  match  the  arrival  times  of  both  first,  second  and  third  arrivals.  Overall,  given  the 
simplicity  of  our  models,  the  complexity  of  the  crustal  structure  in  this  region  and  the  high 
frequency  character  of  the  data,  the  fit  is  better  than  expected.  Our  split  models  give  a  much  better 
fit  than  the  single  crustal  structure  model,  as  shown  in  Figure  3. 

Discussion  of  Receiver  Models 

An  examination  of  regional  refraction  data  and  teleseismic  arrival  time  data  (Ryaboy, 
Suteau-Henson  1989)  indicates  that  considerable  lateral  heterogeneity  exists  in  both  the  crust  and 
upper  mantle  structure  in  this  region,  corresponding  to  the  differing  tectonic  provinces  in  which  the 
waves  travel;  in  the  region  north  of  the  test  site,  the  waves  travel  through  a  shield  region;  to  the 
west  they  travel  through  the  Caucasus  region,  while  to  the  south  and  east  they  travel  through  the 
thick  crustal  region  of  the  Tien  Shan.  This  lateral  heterogeneity  shows  up  in  the  various  crustal 
and  upper  mantle  models  for  the  receiver  structures. 

These  differing  crustal  and  upper  mantle  structures  give  significantly  different  arrival  times 
for  the  various  regions.  This  heterogeneity  is  incorporated  in  our  starting  models.  For  the  stations 
NVS  and  FRU,  for  which  no  previous  source  structure  model  exists,  we  use  the  average  of  the 
source  model  and  the  nearest  known  receiver  structure  as  the  starting  model.  For  the  starting 
model  for  station  BRV,  we  use  the  crust  and  upper  mantle  structure  of  Oreshin,  et  al.  (1992)  as 
obtained  by  teleseismic  modeling  of  P  wave  data.  At  station  WMQ,  we  utilize  the  starting  model  of 
Mangino  and  Ebel  (1992)  to  obtain  a  starting  receiver  crustal  structure.  We  then  systematically 
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vary  these  models  to  get  the  best  fitting  velocity  model  in  each  region.  Our  receiver  models  have  a 
fair  amount  of  nonuniqueness;  it  is  difficult  to  determine  exact  receiver  structures  solely  on  the 
basis  of  our  work.  Our  models  represent  the  best  fit  of  modified  preexisting  receiver  models  to  the 
data  at  BRV  and  WMQ,  and  the  best  fit  of  regional  known  crustal  and  upper  mantle  structures  to 
the  data  at  NVS  and  FRU.  We  have  much  more  confidence  in  the  accuracy  of  our  BRV  and  WMQ 
models,  for  which  we  have  both  quality  data  and  accurate  starting  models,  than  at  NVS  and  FRU, 
for  which  we  do  not. 

Our  final  models  show  considerable  similarity;  however,  the  location  of  the  crust  mantle 
boundary  varies  from  location  to  location,  as  does  the  upper  mantle  transition  depth  from  8.30 
km/sec  to  8.55  knVsec.  At  station  NVS,  the  crust  has  a  thickness  of  50  km  depth  with  a  lower 
crustal  layer  of  7.20  km/sec  lower  crustal  layer  between  30  and  50  km  depth.  The  upper  mantle 
velocity  in  the  region  between  50  and  70  km  has  a  velocity  of  8.23  km/sec,  a  velocity  of  8.35 
km/sec  in  the  region  between  70  and  90  km  depth,  and  a  velocity  of  8.50  km/sec  in  the  region 
between  90  and  1 10  km  depth. 

At  station  FRU,  the  Pn  first  arrival  indicates  a  high  velocity  upper  crustal  layer  of  about  6.0 
km/sec  in  the  upper  10  km  with  a  6.5  km/sec  lower  crustal  layer  of  20  km  thickness  between  10 
and  30  km  depth  and  a  lower  crustal  velocity  of  7.0  km/sec  in  the  region  between  30  and  45 
km.depth.  The  upper  mantle  in  the  region  between  45  and  65  km  depth  has  a  velocity  of  8.25 
km/sec,  a  velocity  of  8.35  km/sec  in  the  region  between  70  and  90  km  depth,  a  velocity  of  8.55 
km/sec  in  the  region  between  90  and  1 10  km  depth,  and  a  velocity  of  8.65  km/sec  in  the  region 
between  1 10  and  250  km  depth. 

At  station  BRV,  our  model  is  similar  to  that  of  Qreshin  et  al.  ( 1992).  We  find  that  the  crust 
has  a  thickness  of  about  50  km  with  a  mid-crustal  velocity  of  6.5  km/sec  and  a  lower  crustal 
velocity  of  6.8  km/sec,  with  a  sharp  transition  at  50  km  depth  to  an  upper  mantle  velocity  of  8.25 
km/sec  extending  down  to  about  90  km.  In  the  region  between  90  and  1 10  km  depth,  the  mantle 
has  a  velocity  of  8.37  km/sec,  and  a  velocity  of  8.25  km/sec  between  1 10  and  130  km  depth. 
Below  130  km  depth  the  mantle  has  a  velocity  of  8.60  km/sec. 

At  station  WMQ,  our  model  is  similar  to  that  of  Mangino  and  Ebel.  The  crustal  depth 
under  station  WMQ  is  55  km,  with  a  low  velocity  crustal  gradient  in  the  upper  10  km,  a  mid  crustal 
velocity  of  6.5  km/sec  in  the  region  between  10  and  30  km  depth,  and  a  lower  crustal  velocity  of 
7.0  km/sec.  The  upper  mantle  has  a  velocity  of  8.25  km/sec  down  to  75  km  depth,  a  lower  upper 
mantle  velocity  of  8.37  km/sec  down  to  95  km,  a  low  velocity  zone  of  8.25  km/sec  between  95 
and  120  km  depth,  a  velocity  of  8.55  km/sec  down  to  180  km,  and  a  velocity  of  8.70  km/sec 
below  180  km  depth.  The  low  velocity  zone  was  needed  in  order  to  match  the  observed  arrival 
times  in  the  late  phases  of  the  data. 
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In  comparison  to  earlier  work,  we  find  that  most  of  our  crustal  structures  show  the  greatest 
resemblance  to  the  U.S.S.R.  shield  models  of  Riviere-Barbier  et  al.  (1992),  which  show  a  crust 
mantle  boundary  at  about  43  km  depth  in  the  shield  regions  in  the  C.I.S.  and  a  transition  from 
about  8.30  km/sec  to  about  8.55  km/sec  at  90  km  depth.  Our  models  somewhat  resemble  the 
U.S.S.R.  shield  models  of  Burdick  et  al.  (1991);  however,  as  they  are  unable  to  exactly  identify 
the  correct  model  for  each  region,  comparison  with  their  work  is  difficult.  Our  BRV  model  better 
resolves  the  details  of  the  upper  mantle  structures  of  Oreshin  et  al.  (1992),  and  our  WMQ  model 
better  resolves  the  details  of  the  mantle  than  the  structure  of  Mangino  and  Ebel  (1992).  All  of  our 
receiver  structures  incorporate  a  layered  crustal  model  in  the  upper  10  km  to  account  for  the 
gradient  in  the  velocity  in  the  upper  crust,  as  shown  in  previous  source  studies. 

Source  Depth  from  Regional  Seismograms 

One  of  the  questior.s  raised  by  our  research  is  whether  our  split-source  receiver  model 
offers  the  possibility  of  determining  the  depth  of  an  event,  and  whether  our  results  are  useful  in 
solving  the  problem  of  earthquake-explosion  discrimination.  In  order  to  test  the  utility  of  our  split- 
source  receiver  model,  we  computed  synthetic  seismograms  for  2  sets  of  models  in  a  slightly 
simplified  six  layer  crustal  structure  (Figure  6).  We  computed  seismograms  at  a  distance  of  955 
km  for  a  receiver  structure  similar  to  that  of  WMQ  and  at  a  distance  of  685  km  for  a  receiver 
structure  similar  to  that  of  BRV.  We  examine  a  range  of  explosion  and  earthquake  source  depths 
for  each  model.  For  the  explosion  source  model  we  used  source  depths  of  0.2, 0.4,  and  0.8  km 
depth,  typical  of  what  would  be  expected  at  a  drilled  test  site.  For  the  earthquake  source  model, 
we  used  a  source  depth  of  5  km  typical  of  those  found  in  the  upper  crust. 

We  computed  Green's  functions  and  convolved  them  with  two  different  source  time 
functions.  For  the  explosion  source,  we  used  the  standard  Mueller-Muiphy  source  as  reported 
earlier.  For  the  earthquake  source  model,  we  used  a  trapezoidal  time  function  of  about  1  second 
duration  for  a  strike  slip  earthquake  similar  to  that  of  a  magnitude  5  earthquake.  We  low  pass 
filtered  the  resulting  events  with  a  low  pass  filter  with  a  cutoff  of  2  Hz.  We  find  that  for  the 
explosion  model,  there  are  slight  changes  in  the  character  of  the  seismograms  depending  on  source 
depth;  the  deeper  explosion  events  have  a  slightly  greater  high  frequency  component  directly  after 
the  first  P  wave  arrival.  The  earthquake  waveforms  are  not  readily  distinguished  from  the 
explosion  waveforms.  Thus  we  are  uncertain  whether  our  work  can  be  used  definitively  for 
earthquake-explosion  discrimination. 
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Conclusions 

We  computed  generalized  ray  synthetic  seismograms  of  explosions  for  events  recorded  in 
the  C.I.S.  between  1971-1989.  Our  split  source-receiver  model  computes  the  Cagniard  ray 
parameters  using  the  split  source-receiver  phase  arrival  times,  and  the  ray  amplitudes  using  the 
multiplied  reflection  and  refraction  coefficients.  We  were  able  to  produce  seven  and  eight  layer 
models  which  incorporate  regional  refraction  data  and  which  match  both  the  timing  and  relative 
amplitudes  of  the  first  and  second  arrivals  in  the  first  10  seconds  of  the  data  using  our  split  model. 

Our  source  model  was  obtained  from  the  woric  of  Quin  and  Thurber  (1992).  This  model 
had  a  crust-mantle  boundary  at  about  45  km  depth,  an  upper  mantle  velocity  of  8.25  km./sec 
between  45  and  70  km  depth,  a  velocity  of  8.37  km/sec  in  the  region  between  70  and  90  km,  a  low 
velocity  zone  of  8.25  km/sec  in  the  region  ^twee.  '  0  and  1 10  km  depth  and  a  velocity  of  8.65  km 
below  130  km.  Our  receiver  models  show  ariety  of  crust  and  upper  mantle  structures.  Stations 
NVS,  BRV,  and  FRU  have  a  crust-mantle  boundary  of  about  47  km  depth,  while  station  WMQ 
has  a  crustal  depth  of  55  km.  Most  receiver  models  had  an  upper  mantle  velocity  of  8.25  km/sec  in 
the  region  between  50  and  70  km  depth  and  a  mid-upper  mantle  velocity  of  8.65  km/sec  in  the 
region  below  about  1 30  km  depth.  Our  modeling  results  may  be  useful  in  improving  our  ability  to 
discriminate  between  earthquakes  and  explosions,  and  in  determining  the  source  depth  of 
earthquakes. 
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BRV  53.0548  70.2763 

FRU  42.8333  74.6167 

NVS  54.9000  83.3000 

WMQ  43.8210  87.6950 
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Figure  Captions 


Figure  1.  Epicenter  and  station  location  map. 

Figure  2.  Illustration  of  ray  paths  for  a  split  source  receiver  model  showing  the  paths  Pn,  PmP, 
Pg,  2 PmP  and  2Pn,  and  other  crustal  and  upper  mantle  phases. 

Figure  3.  Comparison  of  synthetic  seismograms  at  stations  NVS  and  WMQ  computed  using 
single  model  method  and  double  model  method  compared  with  data.  Traces  are,  in  order,  split 
model,  source  only  model,  receiver  only  model,  and  data. 

Figure  4.  Data  recorded  at  stations  FRU,  NVS,  BRV,  and  WMQ  (bottom  in  each  panel)  compared 
with  synthetic  seismograms  (top  in  each  panel). 

Figure  5.  Velocity  models  showing  source  model  and  receiver  models  for  each  station. 

Figure  6.  Synthetic  seismograms  for  explosions  computed  at  depths  of  0.2, 0.4  and  0.8  km  depth 
compared  against  a  strike  slip  earthquake  computed  at  a  depths  of  5  km  (a)  at  a  distance  of  955  km 
for  a  receiver  structure  similar  to  that  of  WMQ  and  (b)  at  a  distance  of  685  km  for  a  receiver 
structure  similar  to  that  of  BRV. 
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Figure  1 
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Ray  Geometry  ol  Split  Source  Receiver  Model 


Figure  2 
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Figure  5 
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Abstract.  Locations  of  20  nuclear  explosions  from  1987 
to  1989  at  the  Balapan  test  site,  Kazakhstan,  are  derived  with  a 
precision  of  about  100  m  from  time-sequence  SPOT  satellite 
images  combined  with  teleseismic  epicenter  estimates.  For 
most  events,  there  is  little  or  no  ambiguity  in  the  association 
between  seismic  events  and  image  features.  The  locations 
determined  by  JED  are  associated  with  formal  error  estimates 
(95%  confidence  ellipses)  that  are  significantly  too  small  In 
16  out  of  19  cases,  the  events  are  found  from  satellite  imagery 
to  lie  outside  these  ellipses.  Possible  causes  are  nonuniform 
observation,  poor  azimuthal  coverage,  and  source  region 
heterogeneity.  Our  results  should  prove  valuable  for  studies 
of  seismic  waveforms,  arrival  time  and  amplitude  tomography, 
event  location,  tectonic  release,  and  seismic  coupling  at 
Balapan,  and  demonstrate  the  utility  of  combined  seismic  and 
satellite  image  analyses  for  seismic  verification  studies. 

Introduction 

The  former  Soviet  Union  conducted  its  last  nuclear 
explosion  at  the  Kazakhstan  (Semipalatinsk)  Test  Site  (KTS), 
in  East  Kazakhstan,  in  October,  1989.  In  August  1991,  the 
newly  independent  country  of  Kazakhstan  closed  the  lest  site 
and  in  May  1992  undertook  to  enter  the  Non-Proliferation  - 
Treaty  as  a  non-nuclear  weapon  state.  Prior  to  its  closing, 
several  hundred  nuclear  explosions  had  been  conducted 
underground  in  Kazakhstan,  making  KTS  second  in  numbers 
of  explosions  only  to  the  Nevada  Test  Site  of  the  US. 

The  signals  of  nuclear  explosions  are  of  great  interest  since 
they  permit  detailed  studies  of  earth  structure  and  properties  in 
the  source  region.  To  address  these  issues,  it  is  often 
necessary  to  obtain  very  accurate  locations  for  the  explosions. 
Perhaps  the  greatest  potential  is  in  "source  array"  studies  (e.g., 
Goldstein  et  al.,  1992;  McLaughlin  et  al.,  1992),  wherein  a 
precise  knowledge  of  source  locations  allows  a  controlled 
study  of  propagation  effects.  The  study  of  tectonic  release 
could  be  advanced  greatly  by  having  accurate  locations  of 
explosions  with  respect  to  nearby  faults.  Many  other  travel 
time,  amplitude,  waveform,  and  source  location  issues  can 
also  be  executed  at  levels  of  precision  not  possible  with 
earthquake  sources,  once  accurate  explosion  locations  are 
available.  In  this  paper,  we  give  esdmates,  to  a  new  level  of 
precision,  of  the  locations  of  20  nuclear  explosions  in  the 
Balapan  area  of  KTS  between  1987  and  1989. 


Information  on  the  locations  and  characteristics  of  96  KTS 
explosions  prior  to  1973  has  been  made  publicly  available  by 
Bocharov  et  al.  (1989)  and  reported  in  the  U.S.  literature  by 
Vergino  (1989a, b).  Of  these  events,  7  were  from  Balapan 
(Shagan  River)  subregion  of  KTS.  There  have  been  a  number 
of  teleseismic  joint  location  studies,  such  as  Marshall  et  al. 
(1984)  and  Lilwall  and  Farthing  (1990).  The  Marshall  et  al. 
(1984)  study  used  the  location  of  the  January  1965  cratering 
explosion  in  a  LANDSAT  image  to  provide  a  master  event  for 
their  analysis.  Lilwall  and  Farthing  (1990)  included  nearly  all 
of  the  events  listed  by  Bocharov  et  al.  (1989)  as  master  events. 

We  report  the  results  of  a  simple  and  direct  approach  to  the 
determination  of  accurate  locations  for  recent  (post  1986)  KTS 
explosions  in  Balapan  by  combining  time-sequence  satellite 
images  of  KTS  with  the  known  location  information  and 
teleseismic  location  estimates.  High-resolution  images  of 
Balapan  from  Satellite  Pour  l'Observation  de  la  Terre  (SPOT) 
are  available  at  irregular  intervals  between  1986  and  late  1989: 
860617,  870807  (partial  coverage),  880826, 881011. 890603, 
890708  (partial  coverage),  and  89101.  Fresh  features  can  be 
associated  with  explosions  occurring  in  the  corresponding  lime 
period.  We  estimate  our  final  location  errors  to  be  about  100 
m,  based  on  comparisons  between  our  locations  and  (a)  those 
available  for  a  small  number  of  events,  and  (b)  those  picked 
independently  by  one  of  us  from  two  SPOT  images  available 
at  the  DARPA  Center  for  Seismic  Studies  in  Arlington,  VA. 
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Fig.  1.  JED  locaiions  of  the  most  recent  20  Balapan 
nuclear  explosions  (post  1986)  with  the  19 88  JVE  used  as 
the  master  event.  The  dashed  boxes  indicate  the  areas  of 
the  SPOT  images  shown  in  Fig  2. 
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Tahle  1.  JED  locations  for  Balapan  explosions.  1987-  1989 
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JED  Analyses  of  Balapan  Explosions 

The  studies  of  Marshall  et  al.  (1984)  and  Lilwall  and 
Farthing  (1990)  both  utilized  the  JED  method  of  Douglas 
(1967)  to  determine  locations  of  KTS  explosions  using  one  or 
more  master  events.  The  Marshall  et  al.  (1984)  study  included 
61  explosions  at  Balapan  through  19X2,  and  used  just  the 
1965  cratering  event  (event  650 1 1 5)  as  a  master  event.  The 
Lilwall  and  Farthing  (1990)  study  included  an  additional  40 
Balapan  events  from  1983  through  the  cessation  of  testing  in 
1989,  and  used  all  7  Balapan  events  listed  by  Bocharov  et  al. 
(1989)  as  master  events.  The  consistency  between  the 
locations  of  the  61  common  events  in  the  two  studies  is  quite 
good,  with  an  RMS  difference  of  only  about  1  km.  Or  own 
analysis  using  both  JED  and  algorithm  MLOC  (E.  Bergman, 
personal  comm.,  1992)  indicate  the  locations  are  stable  within 
about  1  km.  The  main  point  of  this  paper  is  that  with  satellite 
images  we  can  improve  this  precision  by  about  a  factor  of  10. 

Preliminary  locations  of  the  20  post-1986  explosions  at 
Balapan  and  their  estimated  uncertainties  are  shown  in  Figure 
1  and  Table  1,  computed  using  algorithm  JED  (Douglas, 

1967)  and  ISC  data,  with  the  epicenter  of  the  JVE  fixed  at  its 
image  location.  Event  depths  were  fixed  by  scaling  to 
magnitude  with  the  formula  of  Jih  and  Wagner  (1991), 
log(DOB)  =  0.31  mb  +  0.835  (Table  1).  These  locations 
provide  the  starting  point  for  analysis  of  the  SPOT  images. 

Satellite  Image  Analysis 

Nine  events  from  the  Bocharov  et  al.  (1989)  list  were 
employed  as  control  points  for  image  rectification,  7  from 
Balapan  and  2  nearby  throwout  craters  known  as  T- 1  and  T-2 
(Leith  and  Simpson,  1990).  However,  we  suspect  that  the 
location  of  event  7  is  in  error  in  Bocharov  et  al.  (1989),  as 
there  is  a  cratered  feature  at  that  location  (possibly  event  10, 
mb  =  4.4, 740416).  Murphy  and  Jenab  (1992)  associate  event 


7  with  the  same  faint  feature  chosen  by  us  (at  50.0356, 
79.0108),  about  1.5  km  northeast  of  the  location  provided  by 
Bocharov  et  al.  (1989).  The  rectification  produced  an 
excellent  fit  to  the  control  points  (RMS  misfit  of  about  50  m). 

Identifying  the  sites  of  the  post- 1986  explosions  involved 
the  comparison  of  the  time-sequence  SPOT  images  with  the 
estimated  event  locations  indicated  in  Table  1.  Figure  2  shows 
"before-and-after"  images  acquired  on  860617  and  891015  of 
areas  containing  identified  explosion  sites.  The  images  are  8 
km  by  6  km,  and  their  locations  are  indicated  in  Figure  1.  The 
coordinates  and  event  associations  are  presented  in  Table  2. 

Figures  2a  and  b  show  the  1988  JVE  (event  94,  880914) 
and  the  adjacent  events  86  (870802)  and  96  (881217).  SPOT 
images  exist  both  before  and  after  each  of  these  events,  so 
these  associations  are  quite  certain.  Coordinates  for  the  JVE 
have  been  published  by  Murphy  and  Jenab  (1992);  their 
location  of  49.8788®  N,  78.8225°  E  agrees  well  with  our 
estimate  of  49.8781°  N,  78,8239°  E  (less  than  130  m  apart). 
Figures  2  c  and  d  cover  the  area  just  north  of  the  JVE,  with 
events  90  (880213),  91  (880403),  and  101  (891019).  Event 
101  post-dates  the  last  SPOT  image,  but  the  preparation  area 
for  the  explosion  (on  October  19)  is  visible  in  the  image  of 
891015.  Because  events  90  and  91  are  well-separated  both 
from  each  other  and  from  other  events  close  in  time,  it  is 
unlikely  that  their  identities  are  in  error.  Figures  2  e  and  f 
cover  far  NE  Balapan  with  events  93  (880614),  95  (881 1 12), 
and  100(890902).  There  is  complete  before-and-after 
coverage  for  these  events.  The  difficulty  of  distinguishing  the 
emplacement  point  for  100  from  a  probable  nearby  trailer  area 
lend  some  uncertainty  to  this  association.  Figures  2  g  and  h 
cover  the  area  just  west  of  the  JVE  with  events  84  (870417), 

87  (871115),  89  (871227),  and  99  (890708).  The  lack  of 
coverage  of  the  western  area  in  the  800807  image  hinders  our 
association  effort.  The  feature  identified  as  shot  84  is  present 
in  the  860617  image,  though  its  appearance  has  changed  in  the 
next  available  image  (880826).  Events  87  and  89  are  far 
enough  apart  from  each  other  and  other  nearby  events  that  their 
identification  is  reasonably  certain  despite  their  closeness  in 
time.  For  event  99,  we  have  excellent  before-and-after 
coverage,  so  its  identification  is  quite  certain.  Our  location  for 
this  event  of  49.8675  N,  78.7792  E  is  only  about  70  m 
different  from  the  location  of  49.8681  N,  78.7792  E  provided 
to  us  by  Adushkin  et  al.  (1992)  after  we  had  prepared  our 
Table  2.  Figures  2  i  and  j  cover  the  area  just  north  of  the  JVE 
with  events  82  (870312)),  88  (871213),  and  97  (890122). 
Again,  because  event  82  is  spatially  well-separated  from  other 
events  close  in  time,  it  is  unlikely  that  its  identity  is  in  error. 
We  have  complete  before-and-after  coverage  for  events  88  and 
97,  so  these  associations  are  quite  certain.  Figures  2  k  and  1 
cover  the  area  NW  of  the  JVE  with  events  83  (870403),  85 
(870620),  92  (880504),  and  98  (890212).  Our  associations 
for  events  83  and  98  have  high  confidence;  note  also  the 
apparent  observation  point  to  the  ENE  of  event  98.  For  event 
85,  it  is  possible  that  the  shot  point  is  actually  the  smaller 
feature  ESE  of  our  pick;  we  interpret  that  to  be  an  observation 
point  Our  association  for  event  92  is  somewhat  less  certain, 
but  there  is  a  relatively  fresh  feature  visible  in  multi-spectral 
images  that  is  not  readily  apparent  in  panchromatic  images. 

Discussion  and  Conclusions 

Event  origin  times  have  been  computed  relative  to  that  of 
event  99  (03  h  46  m  57.64  s)  as  given  by  Adushkin  et  al. 
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Fig.  2.  Beforc-and-after  SPOT  images  for  explosion  sites.  Each  pair  shows  images  from  860617  (top),  and 
891015  (bottom)  and  the  explosion  site,  as  indicated  below.  a,b)  The  area  of  events  880914  (the  JVE;  above  left), 
870802  (below  center),  and  881217  (above  right),  c.d)  The  area  just  NE  of  the  JVE,  with  events  880213  (below 
center),  880403  (below  center),  and  891019  (below  center).  e,f)  The  far  NE  part  of  Balapan,  with  events  880614 
(below  center),  881112  (above  center),  and  890902  (above  center),  g.h)  The  area  just  W  of  the  JVE,  with  events 
870417  (below  left),  871115  (above  right),  871227  (above  center),  and  890708  (above  center).  i,j)  The  area  just  N 
of  the  JVE,  with  events  870312  (above  center),  871213  (above  center),  and  890122  (below  center).  k,l)  The  area 
NW  of  the  JVE,  with  events  870403  (above  center),  870620  (above  center),  880504  (below  center),  and  890212 
(above  center).  ©  1993  CNES.  Provided  by  SPOT  Image  Corporation. 


(1992),  keeping  all  epicenters  and  depths  fixed  at  the  values 
given  in  Table  2.  The  resulting  values  (Table  2)  average  2.4  s 
earlier  than  those  estimated  by  Lilwall  and  Farthing  (1990). 

In  Figure  3,  the  locations  derived  from  satellite  image 
analysis  (Table  2  and  Figure  2  a-1)  are  compared  to  the  JED 
locations  (Table  1  and  Figure  1).  The  mislocation  vectors 
trend  NNE  in  NE  Balapan  versus  SSW  in  SW  Balapan, 
averaging  1.2  km  for  the  19  unconstrained  events,  compared 
to  an  average  error  ellipse  dimension  of  0.7  km.  The  location 
derived  from  SPOT  images  falls  within  the  95%  confidence 
ellipse  for  only  3  events.  Doubling  the  size  of  the  error 
ellipses  increases  this  to  9  events.  Clearly,  JED  seriously 
underestimates  event  mislocation.  Nonuniform  observation 


(Pavlis,  1992),  poor  azimuthal  coverage,  and  near-source 
velocity  heterogeneity  are  possible  causes  of  this  discrepancy. 
Each  could  produce  systematic  errors  that  are  not  accounted 
for  in  the  JED  estimation  of  uncertainty,  which  assumes  that 
errors  are  uncorrelated.  These  hypotheses  could  be  tested 
once  a  larger  number  of  accurate  locations  are  determined 
Our  final  event  locations  have  value  for  a  variety  of  studies. 
Methods  for  improving  location  and  uncertainty  estimates 
could  be  tested.  They  can  be  used  as  controlled  sources  for 
waveform  modeling  studies  (Goldstein  et  a!.,  1992),  or  to 
improve  tomographic  studies  of  KTS  (McLaughlin  et  al. 
(1992).  Characteristics  of  the  emplacement  points  might  ■ 
provide  information  on  yield  and  depth  of  burial.  The  relation 
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Table  2.  Locations  of  Balanan  explosions  from  SPOT  images 
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09  49  57.4 

49.9300 

2.8,94,56 

♦JVE 


Fig.  3.  Comparison  of  SPOT  image  (squares)  and  JED 
teleseismic  (diamonds)  locations  for  the  recent  20  nuclear 
explosions  at  Balapan.  The  JVE  was  used  as  the  master 
event.  In  16  out  of  19  cases,  mislocations  exceed  the  95% 
confidence  bounds  indicated  in  Fig.  1. 


between  explosion  locations  and  regional  faults  or  preexisting 
rubble  zones  could  provide  important  information  for  studies 
of  tectonic  release  or  seismic  coupling.  These  types  of 
analyses  should  prove  useful  for  the  monitoring  of  nuclear 
explosions  at  other  historic  and  possibly  future  test  sites. 
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2.4.4.  Accurate  locations  of  nuclear  explosions  at  Balapan, 

Kazakhstan,  1973  to  1985 

Introduction 

Over  the  years  1965  through  1989,  the  former  Soviet  Union  is  thought  to  have  conducted 
101  nuclear  tests  at  the  Shagan  River  test  site  (Lilwall  and  Farthing,  1990).  All  of  these  events 
were  recorded  and  located  teleseismically  by  a  large  number  of  stations.  Exact  location  of  these 
events  using  teleseismic  means  alone  is  impossible  due  to  the  uncertainties  inherent  in  seismic 
location  techniques.  However,  in  the  years  1974  to  1982  the  Soviet  test  site  was  imaged  often  by 
the  LANDSAT  satellite,  and  imaged  frequently  in  the  years  between  1986-1989  by  the  SPOT 
satellite.  These  satellite  images  show  clearly  identifiable  shot  areas;  in  most  cases,  a  circular  area  is 
visible  on  the  satellite  image  at  or  near  the  junction  of  two  or  more  roads.  The  identification  of  shot 
areas  by  previous  work  on  SPOT  imagery  and  teleseismic  estimates  of  shot  location  has  proven  to 
be  an  effective  method  of  obtaining  exact  locations  of  nuclear  explosions  for  Shagan  River 
(Thurber  et  al.,  1993). 

In  this  report,  we  present  locations  for  suspected  explosions  from  the  Shagan  River  test  site 
in  the  years  1965  to  1985.  Starting  estimates  of  the  event  locations  were  obtained  from  a  revision 
of  the  results  Lilwall  and  Farthing  (1990),  using  the  program  JED,  in  which  27  known  events 
were  used  to  constrain  the  teleseismic  locations,  compared  to  only  7  used  by  Lilwall  and  Farthing 
(1990).  The  7  master  events  used  by  Lilwall  and  Farthing  (1990)  are  those  events  from  the  years 
1965  to  1972  published  by  Bocharov  et  al.  (1989).  The  additional  20  master  events  we  use  are 
those  whose  locations  were  determined  from  SPOT  images  by  Thurber  et  al.  (1993). 

Using  the  revised  location  estimates,  we  endeavor  to  make  accurate  associations  between 
event  epicenters  and  identifiable  shot  points  on  rectified  LANDSAT  and  SPOT  images.  The  4- 
band  LANDSAT  images  are  enhanced  by  various  techniques.  This  enhancement  allows  clear 
identification  of  a  large  number  of  new  shot  areas  occurring  in  a  sequence  of  images.  In  this 
manner,  stepping  forward  in  time  from  image  to  image,  overlaying  shots  from  previous  images  on 
top  of  the  current  image,  we  can  identify  virtually  all  explosions  occurring  in  the  years  between 
1973  and  1985. 


LANDSAT  and  SPOT  Data 

The  images  used  in  this  study  consist  of  a  set  of  LANDSAT  images  recorded  in  the  years 
1974  to  1982.  Examples  are  shown  in  Figure  7.  These  data  are  recorded  at  a  spatial  resolution  of 
80  meters,  which  was  sufficient  to  resolve  the  location  of  the  shot  areas.  The  images  were 
registered  to  a  common  set  of  ground  control  points  obtained  from  a  rectified  SPOT  1986  image. 
We  were  able  to  register  each  of  the  images  to  within  about  one  to  two  hundred  meters  of  each 
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other.  Each  registered  image  covered  the  area  between  78.6°  and  79.1°  E  in  longitude  and  49.7° 
and  50. 1°  N  in  latitude.  This  covered  the  estimated  teleseismic  location  of  all  explosions  with  a 
magnitude  of  about  5.0  and  above  at  the  Balapan  test  site. 

This  common  rectification  enabled  us  to  make  accurate  comparisons  of  shot  location 
between  various  images.  By  overlying  the  identified  shot  points  from  one  image  on  another,  we 
can  make  accurate  comparisons  of  shot  points  from  one  year  to  the  next.  The  emplacement  point 
identifications  were  made  on  the  rectified  images  by  enhancing  the  images  to  increase  the  contrast 
of  the  shot  points  from  the  surrounding  regions.  On  the  enhanced  images,  we  can  make  accurate 
picks  on  which  shot  points  occurred  since  the  last  image  was  acquired.  In  this  manner,  proceeding 
forward  in  time  image  by  image  since  1974,  by  matching  new  emplacement  locations  with 
available  teleseismic  locations  we  can  accurately  identify  essentially  all  the  shot  emplacement 
points.  Our  final  locations  are  provided  in  Table  3,  including  all  101  events  for  convenience. 

Discussion 

We  find  that  we  are  able  to  associate  teleseismic  short  locations  with  satellite  images  of  shot 
points  for  virtually  all  of  the  shots  occurring  between  1973  and  1985.  For  most  events  the 
teleseismic  locations  and  LANDSAT  shot  point  estimates  agree  to  within  about  1  to  2  km.  In  the 
process  of  carrying  out  this  work,  we  have  learned  that  one  of  the  events  (#72  15-Sep-84)  is 
reported  by  Russian  scientists  to  be  a  large  chemical  explosion,  not  a  nuclear  explosion.  This 
event  is  one  of  the  few  for  which  our  association  is  uncertain.  Others  which  we  deem  somewhat 
uncertain  are  indicated  by  an  asterisk  (*)  in  Table  3;  these  include  events  10, 12,  and  65.  Event  10 
is  one  of  the  smallest  explosions  (mb  =  4.35)  with  among  the  largest  location  uncertainties.  We 
have  associated  it  with  an  apparent  cratering  feature  that  in  fact  is  close  to  (but  not  exactly  equal  to) 
the  location  given  by  Bocharov  et  al.  (1990)  for  event  7.  Thus  this  association  is  probably  to  be 
considered  controversial.  In  the  cases  of  events  12  and  72,  the  regions  in  which  they  are  thought 
to  be  located  (based  on  the  teleseismic  location)  are  quite  bright  and  complex,  so  it  is  difficult  to 
detect  new  features,  especially  at  the  low  sptial  resolution  of  the  LANDSAT  data.  For  event  65,  it 
is  quite  possible  that  the  bright  feature  immediately  southwest  of  our  chosen  point  is  the  explosion 
site,  but  the  latter  feature  appears  to  be  just  a  road  junction,  not  an  explosion  site.  Overall,  the 
combined  use  of  satellite  imagery  with  teleseismic  (and  regional)  locations  should  be  quite  useful  in 
monitoring  a  global  test  ban  for  the  case  in  which  suspected  test  areas  involve  vertical  drilling. 
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3. 


A  plane-layered  velocity  model  has  proven  adequate  for  waveform  modeling  and  event 
location  a  t  near-regional  (<  300  km)  distances  in  the  vicinity  of  the  Kazakhstan  Test  Site  (KTS). 
We  have  verified  that  the  approach  of  Thurber  et  al.  (1989)  and  Li  and  Thurber  (1991)  of  using 
multiple  secondary  body  wave  arrivals  for  event  location  is  justified  by  successfully  identifying 
and  modeling  Pg  and  PmP  secondary  phases,  and  relocating  events  using  these  confirmed  phases. 
This  approach  would  be  practical  in  other  areas  with  some  calibration  shot  or  other  master  event 
information  and  sparse  near-regional  station  coverage. 

For  moderate  regional  distances  (600  to  900  km),  we  have  applied  a  novel  modification  of 
the  generalized  ray  method  to  synthesis  seismograms  in  a  manner  that  accounts  for  differences  in 
structure  between  the  source  region  and  receiver  region.  Models  for  structure  beneath  four  Central 
Asia  stations  (BRV,  FRU,  NVS,  and  WMQ)  have  been  derived  using  events  from  KTS  and  the 
source  region  structure  derived  above. 

We  give  estimates,  to  a  new  level  of  precision,  of  the  locations  of  20  nuclear  explosions  at 
the  Balapan  Test  Site  between  1987  and  1989.  We  have  very  high  confidence  in  the  identification 
from  SPOT  images  of  the  explosion  emplacement  points  for  16  of  the  20  explosions  with  mb  >  5 
post-dating  1986.  Improved  temporal  image  coverage  and/or  availability  of  additional  multispectral 
images  from  other  sources  might  help  resolve  remaining  uncertainties  in  event  identification. 
Comparing  the  site  locations  from  rectified  satellite  images  with  teleseismic  joint  epicenter 
determination  (JED)  results  leads  us  to  conclude  that  seismic  epicenter  uncertainties  are  generally 
underestimated  by  a  factor  of  2  or  more  for  Shagan  River.  We  suggest  that  non-i;  liform 
observation  (Pavlis,  1992)  combined  with  poor  azimuthal  coverage  may  be  a  major  cause  of  this 
discrepancy.  For  the  events  in  the  time  period  1973  to  1985,  we  have  been  successful  in 
identifying  explosion  sites  for  virtually  all  of  teleseimically-detected  events,  using  LANDSAT  and 
SPOT  images. 

Our  final  event  locations  have  potential  value  for  a  variety  of  studies.  Methods  for 
improving  location  estimates  and  uncertainties  (VanDecar  and  Cnosson,  1990;  Pavlis,  1992)  could 
be  tested  with  these  events.  They  could  be  used  as  controlled  sources  for  waveform  modeling 
studies  (Goldstein  et  al.,  1992)  and  tomographic  studies  (McLaughlin  et  al.,  1992).  Perhaps  the 
mb  bias  variations  across  the  Shagan  River  area  (Ringdal  et  al.,  1992)  could  be  accounted  for  in 
part  by  lateral  variations  in  velocity  structure.  The  characteristics  of  the  emplacement  points  might 
provide  some  information  on  explosion  yield  and  depth  of  burial.  These  types  of  analyses  could 
prove  useful  for  the  monitoring  of  explosions  at  other  historic  and  possibly  future  test  sites. 
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Figure  7.  Examples  of  LANDSAT  images  used  for  identifying  explosion  sites,  from  (a)  1974, 
(b)  1977,  (c)  1979,  (d)  1981,  and  (e)  a  SPOT  image  from  1986  [©  1992  CNES.  Provided  by 
SPOT  Image  Corporation.]. 
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